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ABSTRACT 

A theoretical model for quasi-spherical subsonic accretion onto slowly rotating mag- 
netized neutron stars is constructed. In this model the accreting matter subsonically 
settles down onto the rotating magnetosphere forming an extended quasi-static shell. 
This shell mediates the angular momentum removal from the rotating neutron star 
magnetosphere during spin-down episodes by large-scale convective motions. The 
accretion rate through the shell is determined by the ability of the plasma to enter the 
magnetosphere. The settling regime of accretion can be realized for moderate accre- 
tion rates M < M, ^ 4 x 10'^ g/s. At higher accretion rates a free-fall gap above the 
neutron star magnetosphere appears due to rapid Compton cooling, and accretion be- 
comes highly non-stationary. From observations of the spin-up/spin-down rates (the 
angular rotation frequency derivative of , and dcf jdM near the torque reversal) of 
X-ray pulsars with known orbital periods, it is possible to determine the main di- 
mensionless parameters of the model, as well as to estimate the magnetic field of 
the neutron star We illustrate the model by determining these parameters for three 
wind-fed X-ray pulsars GX 301-2, Vela X-1, and GX 1-1-4. The model explains both 
the spin-up/spin-down of the pulsar frequency on large time-scales and the irregular 
short-term frequency fluctuations, which can correlate or anti-correlate with the X- 
ray flux fluctuations in different systems. It is shown that in real pulsars an almost 
iso-angular-momentum rotation law with lo ~ \ /R^, due to strongly anisotropic ra- 
dial turbulent motions sustained by large-scale convection, is preferred. 
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1 INTRODUCTION 

X-ray pulsars are highly magnetized neutron stars in binary systems, accreting matter from a 
companion star. The companion may be a low-mass star overfilling its Roche lobe in which case 
an accretion disc is formed. In the case of a high-mass companion, the neutron star may also 
accrete from the strong stellar wind and depending on the conditions a disc may be formed or 
accretion may take place quasi-spherically. The strong magnetic field of the neutron star disrupts 
the accretion flow at some distance from the neutron star surface and forces the accreted matter to 
funnel down on the polar caps of the neutron star creating hot spots that, if misaligned with the 
rotational axis, make the neutron star pulsate in X-rays. Most accreting pulsars show stochastic 
variations in their spin frequencies as well as in their luminosities. Many sources also exhibit long- 
term trends in their spin-behaviour with the period more or less steadily increasing or decreasing 
and in some sources spin-reversals have been observed. (For a thorough review, see e.g. Bildsten 
1997 and references therein.) 

The best-studied case of accretion is that of thin disc accretion (Shakura & Sunyaev 1973). 
Here the spin-up/spin-down mechanisms are rather well understood. For disc accretion the spin- 
up torque is determined by the specific angular momentum at the inner edge of the disc and can 
be written in the form K^u ~ M yfGMRl (Pringle and Rees, 1972). For a pulsar the inner radius 
of the accretion disc is determined by the Alfven radius Ra, Ra ~ M'^^'' , so Ksu ~ M^^^, i.e. for 
disc accretion the spin-up torque is weakly (almost linearly) dependent on the accretion rate (X- 
ray luminosity). In contrast, the spin-down torque for disc accretion in the first approximation is 

independent of M: K^d f^^/Rl, where R^ = (GM/(a»*)^)'^^ is the corotation radius, oj* is the 

neutron star angular frequency and fx is the neutron star's dipole magnetic moment. In fact, accre- 
tion torques in disc accretion are determined by a complicated disc-magnetospheric interaction, 
see, e.g., Ghosh & Lamb (1979), Lovelace et al. (1995) and the discussion in Kluzniak & Rappa- 
port (2007), and correpsondingly can have a more complicated dependence on the mass accretion 
rate and other parameters. 

Measurements of spin-up/spin-down in X-ray pulsars can be used to evaluate a very important 
parameter of the neutron star - its magnetic field. The period of the pulsar is usually close to the 
equilibrium value Pgq, which is determined by the total zero torque applied to the neutron star, 
K = + = 0. So assuming the observed value oj* = InjPeq, the magnetic field of the neutron 
star in disc-accreting X-ray pulsars can be estimated if M is known. 

In the case of quasi-spherical accretion, which can take place in systems where the optical 
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star underfills its Roche lobe and no accretion disc is formed, the situation is more complicated. 
Clearly, the amount and sign of the angular mometum supplied to the neutron star from the cap- 
tured stellar wind are important for spin-up or spin-down. To within a numerical factor (which 
can be positive or negative, see numerical simulations by Fryxell & Taam (1988), Ruffert (1997), 
Ruffert (1999), etc.), the torque applied to the neutron star in this case should be proportional to 
MojeR^, where ojb = InlPs is the binary orbital angular frequency, Rb = 2GM/(Vl + v^^^)^ is 
the gravitational capture (Bondi) radius, is the stellar wind velocity at the neutron star orbital 
distance, and Vgrb is the neutron star orbital velocity. In real high-mass X-ray binaries the orbital 
eccentricity is non-zero, the stellar wind is variable and can be inhomogeneous, etc., so K^u can be 
a complicated function of time. The spin-down torque is even more uncertain, since it is impossi- 
ble to write down a simple equation like -fx^/Rl any more (Rc has no meaning for quasi-spherical 
accretion; for slowly rotating pulsars it is much larger than the Alfven radius where the angular 
momentum transfer from the accreting matter to the magnetosphere actually occurs). For example, 
if one uses for the braking torque -jj^ jRl, the magnetic field in long-period X-ray pulsars turns 
out very high. We think this is a result of underestimating the braking torque. 

The matter captured from the stellar wind can accrete onto the neutron star in different ways. 
Indeed, if the X-ray flux from the accreting neutron star is sufficiently high, the shocked matter 
rapidly cools down due to Compton processes and freely falls down toward the magnetosphere. 
The velocity of motion rapidly becomes supersonic, so a shock is formed above the magneto- 
sphere. This regime was considered, e.g., by Burnard et al. (1983). Depending on the sign of the 
specific angular momentum of falling matter (prograde or retrograde), the neutron star can spin-up 
or spin-down. However, if the X-ray flux at the Bondi radius is below some value, the shocked 
matter remains hot, the radial velocity of the plasma is subsonic, and the settling accretion regime 
sets in. A hot quasi-static shell forms around the magnetosphere (Davies & Pringle, 1981). Due to 
additional energy release (especially near the base of the shell), the temperature gradient across the 
shell becomes superadiabatic, so large-scale convective motions inevitably appear. The convection 
intitiates turbulence, and the motion of a fluid element in the shell becomes quite complicated. If 
the magnetosphere allows plasma entry via instabilities (and subsequent accretion onto the neutron 
star), the actual accretion rate through such a shell is controlled by the magnetosphere (for exam- 
ple, the shell can exist, but the accretion through it can be weak or even absent altogether). So on 
top of the convective motions, the mean radial velocity of matter toward the magnetosphere, a sub- 
sonic settling, appears. This picture of accretion is realized at relatively small X-ray luminositites, 
L;c < 4 X 10^^ erg/s (see below), and is totally different from what was considered in numerical 
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simulations cited above. If the shell is present, its interaction with the rotating magneto sphere can 
lead to spin-up or spin-down of the neutron star, depending on the sign of the angular velocity 
difference between the accreting matter and the magnetospheric boundary. So in the settling ac- 
cretion regime, both spin-up or spin-down of the neutron star is possible, even if the sign of the 
specific angular momentum of captured matter is prograde. The shell here mediates the angular 
momentum transfer to or from the rotating neutron star. 

One can find several models in the literature (see especially lUarionov & Kompaneets 1990 and 
Bisnovatyi-Kogan 1991), from which the expression for the spin-down torque for quasi-spherically 
accreting neutron stars in the form K^^ MR\a>* M^^^ can be derived. Moreover, the ex- 
pression for the Alfven radius Ra in the case of settling accretion is found to have diff'erent depen- 
dence on the mass accretion rate M and neutron star magnetic moment fi, ~ M'^^^^fi^^^^, than the 
standard expression for disc accretion, ~ M'^^'^f/^'^, so the spin-down torque for quasi-spherical 
settling accretion depends on the accretion rate as K^d M^^^^ (see below). 

To stress the difi'erence between quasi-spherical and disc accretion, it is also instructive to 
rewrite the expression for the spin-down torque using the corotation and Alfven radii as Kgd ~ 
-jj^/ ^RlRl ~ -fi^/Rl(Rc/RA?'^ (see more detail below in Section 4). Since the factor (RJRa)^'^ ~ 
(a>K(RA)/oJ*) can be of the order of 10 or higher in real systems, using a braking torque in the form 
IJ^/RI leads to a strong overestimation of the magnetic field. 

The dependence of the braking torque on the accretion rate in the case of quasi- spherical 
settling accretion suggests that variations of the mass accretion rate (and X-ray luminosity) must 
lead to a transition from spin-up (at high accretion rates) to spin-down (at small accretion rates) at 
some critical value of M (or Ra), that diff'ers from source to source. This phenomenon (also known 
as torque reversal) is actually observed in wind-fed pulsars like Vela X-1, GX 301-2 and GX 1+4, 
which we shall consider below in more detail. 

The structure of this paper is as follows. In Section 2, we present an outline of the theory for 
quasi-spherical accretion onto a neutron star magneto sphere. We show that it is possible to con- 
struct a hot envelope around the neutron star through which accretion can take place and act to 
either spin up or spin down the neutron star. In Section 3, we discuss the structure of the inter- 
change instability region which determines whether the plasma can enter the magnetosphere of 
the rotating neutron star. In Section 4 we consider how the spin-up/spin-down torques vary with a 
changing accretion rate. In Section 5, we show how to determine the parameters of quasi- spherical 
accretion and the neutron star magnetic field from observational data. In Section 6, we apply our 
methods to the specific pulsars GX 301-2, Vela X-1 and GX 1-1-4. In Section 7 we discuss our 
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results and, finally, in Section 8 we present our conclusions. A detailed gas-dynamic treatment of 
the problem is presented in four appendices, which are very important to understand the physical 
processes involved. 

2 QUASI-SPHERICAL ACCRETION 

2.1 The subsonic Bondi accretion shell 

We shall here consider the torques applied to a neutron star in the case of quasi-spherical accretion 
from a stellar wind. Wind matter is gravitationally captured by the moving neutron star and a 
bow-shock is formed at a characteristic distance R ~ Rb, where Rb is the Bondi radius. Angular 
momentum can be removed from the neutron star magnetosphere in two ways — either with matter 
expelled from the magnetospheric boundary without accretion (the propeller regime, Illarionov & 
Sunyaev 1975), or via convective motions, which bring away angular momentum in a subsonic 
quasi-static shell around the magnetosphere, with accretion (the settling accretion regime). 

In such a quasi-static shell, the temperature will be high (of the order of the virial temperature, 
see Davies and Pringle (1981)), and the important point is whether hot matter from the shell can 
in fact enter the magnetosphere. Two-dimensional calculations by Eisner and Lamb (1977) have 
shown that hot monoatomic ideal plasma is stable relative to the Rayleigh-Taylor instability at 
the magnetospheric boundary, and plasma cooling is thus needed for accretion to begin. However, 
a closer inspection of the 3 -dimensional calculations by Arons and Lea (1976a) reveals that the 
hot plasma is only marginally stable at the magnetospheric equator (to within 5% accuracy of 
their calculations). Compton cooling and the possible presence of dissipative phenomena (mag- 
netic reconnection etc.) facilitates the plasma entering the magnetosphere. We will show that both 
accretion of matter from a hot envelope and spin-down of the neutron star is indeed possible. 

2.2 The structure of the shell around a neutron star magnetosphere 

To a zeroth approximation, we can neglect both rotation and radial motion (accretion) of matter in 
the shell and consider only its hydrostatic structure. The radial velocity of matter falling through 
the shell Ur is lower than the sound velocity Cj. Under these assumptions, the characteristic cool- 
ing/heating time-scale is much larger than the free-fall time-scale. 

In the general case where both gas pressure and anisotropic turbulent motions are present, Pas- 
cal's law is violated. Then the hydrostatic equilibrium equation can be derived from the equation 
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of motion Eq. (IA16I) with stress tensor components Eq. ( IA19I ) - Eq. (IA21I) and zero viscosity (see 
Appendix A for more detail): 

1 dP, 1 d{P[R^) ^ GM _ 

Here = pc^ly is the gas pressure, and P' stands for the pressure due to turbulent motions: 

P[=p<ul>= pmlc] = yPgni] (2) 

P'^= p <u\>= pm\c] = yPgm\ (3) 

(< >=< uj > +2 < u]_ > is the turbulent velocity dispersion, and are turbulent Mach 
numbers squared in the radial and tangential directions, respectively; for example, in the case of 
isotropic turbulence mj = m\ = {l/3)mj where is the turbulent Mach number). The total 
pressure is the sum of the gas and turbulence terms: Pg + P, = Pg(l + ymj). 

We shall consider, to a first approximation, that the entropy distribution in the shell is constant. 
Integrating the hydrostatic equilibrium equation Eq. ([T]), we readily get 



(y-l\GM 



l^m \ y I R 



1 



y-lGM 

L — ^(7,m,). (4) 

y R 



1 + ym\ - 2{y - l)(mj - ml) ^ 
(In this solution we have neglected the integration constant, which is not important deep inside the 
shell. It is important in the outer part of the shell, but since the outer region close to the bow shock 
at ~ Rb is not spherically symmetric, its structure can be found only numerically). 

Note that taking turbulence into account somewhat decreases the temperature within the shell. 
Most important, however, is that the anisotropy of turbulent motions, caused by convection in the 
stationary case, changes the distribution of the angular velocity in the shell. Below we will show 
that in the case of isotropic turbulence, the angular velocity distribution within the shell is close 
to the quasi-Keplerian one, a){R) ~ R~^^^. In the case of strongly anisotropic turbulence caused 
by convection, mj » , an approximately iso-angular-momentum distribution, oj{R) ~ R~^ is 
realized within the shell. Below we shall see that teh analysis of rela X-ray pulsars favors the 
iso-angular-momentum rotation distribution. 

Now, let us write down how the density varies inside the quasi-static shell for R «; R^. For a 
fully ionized gas with y = 5/3 we find: 
/R \^'^ 

p(R)=p(RA)[-jj (5) 

and for the gas pressure: 

/R \^'^ 

P(R) = P{RA)[-jj . (6) 
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The above equations describe the structure of an ideal static adiabatic shell above the magneto- 
sphere. Of course, at i? ~ i?B the problem is essentially non-spherically symmetric and numerical 
simulations are required. 

Corrections to the adiabatic temperature gradient due to convective energy transport through 
the shell are calculated in Appendix C. 



2.3 The Alfven surface 

At the magnetospheric boundary (the Alfven surface), the total pressure (including isotropic gas 
pressure and the possibly anisotropic turbulent pressure) is balanced by the magnetic pressure 

Pg + P, = Pg(RA)(l + ym?) = . (7) 

The magnetic field at the Alfven radius is determined by the dipole magnetic field and by 
electric currents flowing on the Alfvenic surface 

P^(R^) = —^—-^[^] =^ (8) 
(l+ymf)Sn\RAl Hm 

where the dimensionless coefficient K2 takes into account the contribution from these currents 
and the factor 1/(1 + ym,) is due to the turbulent pressure term. For example, in the model by 
Arons and Lea (1976a, their Eq. 31), K2 = {2.1 Sf w 7.56. At the magnetospheric cusp (where the 
magnetic force line is branched), the radius of the Alfven surface is about 0.51 times that of the 
equatorial radius (Arons and Lea, 1976a). Below we shall assume that Ra is the equatorial radius 
of the magnetosphere, unless stated otherwise. 

Due to the interchange instability, the plasma can enter the neutron star magnetosphere. In the 
stationary regime, let us introduce the accretion rate M onto the neutron star surface. From the 
continuity equation in the shell we find 



P(~Ra) = 7— TTTT^ (9) 



47TUr(Ra)RI 



Clearly, the velocity of absorption of matter by the magnetosphere is smaller than the free-fall 
velocity, so we introduce a dimensionless factor f(u) = ur/ ^2GM/R < 1. Then the density at the 

magnetospheric boundary is 
M 

P(Ra) = , . (10) 

Anfiu) ^1GMIRaR\ 

For example, in the model calculations by Arons & Lea (1976a) f{u) « 0.1; in our case, at high 
X-ray luminosities the value of f{u) can attain » 0.5. It is possible to imagine that the shell is 
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impenetrable and that there is no accretion through it, M ^ 0. In this case ur 0, f{u) 0, while 
the density in the shell remains finite. In some sense, the matter leaks from the magnetosphere 
down to the neutron star, and the leakage can be either small (M 0) or large (M 0). 

Plugging p{R) into Eq. ([8]) and using Eq. dlj) and the definition of the dipole magnetic moment 



(where i?o is the neutron star radius), we find 

= [^z ^<"'^' , T . (11) 

[(r- l)iA(r,m,)(l +ymf)Myf2GM_ 
It should be stressed that in the presence of a hot shell the Alfven radius is determined by 
the static gas pressure at the magneto spheric boundary, which is non-zero even for a zero-mass 
accretion rate through the shell, so the appearance of M in the above formula is strictly formal. 



2.4 Angular momentum transfer 

We now consider a quasi- stationary subsonic shell in which accretion proceeds onto the neutron 
star magnetosphere. We stress that in this regime, i.e. the settling regime, the accretion rate onto 
the neutron star is determined by the denisity at the bottom of the shell (which is directly related 
to the density downstream the bow shock in the gravitational capture region) and the ability of the 
plasma to enter the magnetosphere through the Alfven surface. 

The rotation law in the shell depends on the treatment of the turbulent viscosity (see Appendix 
A for the Prandtl law and isotropic turbulence) and the possible anisotropy of the turbulence due 
to convection (see Appendix B). In the last case the anistropy leads to more powerful radial tur- 
bulence than in the perpendicular directions. Thus, as shown in Appendix A and B, there is a 
set of quasi-power-law solutions for the radial dependence of the angular rotation velocity in a 
convective shell. We shall consider a power-law dependence of the angular velocity on radius, 

co(R) ~ R'" (12) 

We will study in detail the quasi-Keplerian law with n = 3/2, and the iso-angular-momentum 
distribution with n = 2, which in some sense are limiting cases among the possible solutions. 

When approaching the bow shock, R Rb, co ojb- Near the bow shock the problem is not 
spherically symmetric any more since the flow is more complicated (part of the flow bends across 
the shell), and the structure of the flow can be studied only using numerical simulations. In the 
absence of such simulations, we shall assume that the iso-angular-momentum distribution is valid 
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up to the nearest distance to the bow shock from the neutron star which we shall take to be the 
gravitational capture radius Rb, 



where is the stellar wind velocity at the neutron star orbital distance, and Vgrb is the neutron star 
orbital velocity. 

This means that the angular velocity of rotation of matter near the magnetosphere will be 
related toajs via 



(Here the numerical factor to > I takes into account the deviation of the actual rotational law 
from the value obtained by using the assumed power-law dependence near the Alfven surface; see 
Appendix A for more detail.) 

Let the NS magnetosphere rotate with an angular velocity oj* = InjP* where P* is the neutron 
star spin period. The matter at the bottom of the shell rotates with an angular velocity cjm, in general 
different from of. If of > cOm, coupling of the plasma with the magnetosphere ensures transfer of 
angular momentum from the magnetosphere to the shell, or from the shell to the magnetosphere if 
of < (Om- In the general case, the coupling of matter with the magnetosphere can be moderate or 
strong. In the strong coupling regime the toroidal magnetic field component 5, is proportional to 

the poloidal field component Bp as Bp(ojm - oj*)t, and can grow to ~ \Bp\. This regime 

can be expected for rapidly rotating magnetopsheres when oj* is comparable to or even greater 
than the Keplerian angular frequency ojKiRA)', in the latter case the propeller regime sets in. In 
the moderate coupling regime, the plasma can enter the magnetosphere due to instabilities on a 
timescale shorter than that needed for the toroidal field to grow to the value of the poloidal field, 
so B, < Bp. 

2.4.1 The case of strong coupling 

Let us first consider the strong coupling regime. In this regime, powerful large-scale convective 
motions can lead to turbulent magnetic field diffusion accompanied by magnetic field dissipation. 
This process is characterized by the turbulent magnetic field diffusion coefficient 77,. In this case 
the toroidal magnetic field (see e.g. Lovelace et al. 1995 and references therein) is 



RB^2GMI(yl + vU 




(13) 



Bf = — (ojm — )Bp . 



(14) 
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The turbulent magnetic diffusion coefficient is related to the kinematic turbulent viscosity as 77, ^ 
Vf The latter can be written as 

Vt =< Uth > . (15) 

According to the phenomenological Prandtl law which relates the average characteristics of a 
turbulent flow (the velocity Ut, the characteristic scale of turbulence It and the shear ojm - of) 

Ut - lt\C0m- 0J*\. (16) 

In our case, the turbulent scale must be determined by the largest scale of the energy supply to 
turbulence from the rotation of a non-spherical magnetospheric surface. This scale is determined 
by the velocity diff"erence of the solidly rotating magneto sphere and the accreting matter that is 
still not interacting with the magnetosphere, i.e. /, - Ra, which determines the turn-over velocity 
of the largest turbulence eddies. At smaller scales a turbulent cascade develops. Substituting this 
scale into equations Eq. (fT4)) -Eq. (fT6l) above, we find that in the strong coupling regime Bt - Bp. 
The moment of forces due to plasma-magnetosphere interactions is applied to the neutron star and 
causes spin evolution according to: 

loj* = -^mdS = +K{e)K2^ (17) 

where / is the neutron star's moment of inertia, vj is the distance from the rotational axis and K{6) 
is a numerical coefficient depending on the angle between the rotational and magnetic dipole axis. 
The coefficient K2 appears in the above expression for the same reason as in Eq. ([8]). The positive 
sign corresponds to positive flux of angular momentum to the neutron star {co^ > of). The negative 
sign corresponds to negative flux of angular momentum across the magnetosphere (a»„, < of). 

At the Alfven radius, the matter couples with the magnetosphere and acquires the angular 
velocity of the neutron star. It then falls onto the neutron- star surface and returns the angular 
momentum acquired at Ra back to the neutron star via the magnetic field. As a result of this 
process, the neutron star spins up at a rate determined by the expression: 

loj* = +zMR\o)* (18) 

where z is a numerical coefficient which takes into account the angular momentum of the falling 
matter. If all matter falls from the equatorial equator, z = 1 ; if matter falls strictly along the spin 
axis, z = 0. If all matter were to fall across the entire magnetospheric surface, then z = 2/3. 
Ultimately, the total torque applied to the neutron star in the strong coupling regime yields 



,2 



lof = +K{e)K2^ + zMR\cif (19) 
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Using Eq. ([TT]) . we can eliminate M in the above equation to obtain in the spin-up regime 

{OJm > CO*) 

4yf(u) 



. , K(e)K2iu 
lo) = 



4 



l+z- 



(a 



where = GM/{co*)^ is the corotation radius. In the spin-down regime (cOm < to*) we find 



Ico = - 



Ayf{u) 

1 - z- 



(20) 



(21) 



^{y-m+ym^)ifj{y,m,)K{e) 
Note that in both cases must be smaller than Rc, otherwise the propeller effect prohibits ac- 
cretion. In the propeller regime Ra > Rc, matter does not fall onto the neutron star, there are no 
accretion- generated X-rays from the neutron star, the shell rapidly cools down and shrinks and the 
standard lUarionov and Sunyaev propeller (1975), with matter outflow from the magnetosphere is 
established. 

In both accretion regimes (spin-up and spin-down), the neutron star angular velocity co* almost 
approaches the angular velocity of matter at the magneto spheric boundary, co* cOmiRA)- The 
diff'erence between co* and co„, is small when the second term in the square brackets in Eq. (I20l) 
and Eq. (|2T]) is much smaller than unity. Also note that when approaching the propeller regime 
{Ra Rc), the accretion rate decreases, f(u) 0, the second term in the square brackets vanishes, 
and the spin evolution is determined solely by the spin-down term -K{d)f/ /R^^. (In the propeller 
regime, co,„ < cok{Ra), <^m < oj*, co* > cok(Ra) )■ So the neutron star spins down to the Keplerian 
frequency at the Alfven radius. In this regime, the specific angular momentum of the matter that 
flows in and out from the magnetosphere is, of course, conserved. 

Near the equilibrium accretion state (co* ~ cOm), relatively small fluctuations in M across the 
shell would lead to very strong fluctuations in co* since the toroidal field component can change 
its sign by changing from +Bp to -Bp. This property, if realized in nature, could be the distinctive 
feature of the strong coupling regime. It is known (see eg.g. Bildsten et al. 1997, Finger et al. 
201 1) that real X-ray pulsars sometimes exhibit rapid spin-up/spin-down transitions not associated 
with X-ray luminosity changes, which may be evidence for them temporarily entering the strong 
coupling regime. It is not excluded that the triggering of the strong coupling regime may be due to 
the magnetic field frozen into the accreting plasma that has not yet entered the magnetosphere. 



2.4.2 The case of moderate coupling 

The strong coupling regime considered above can be realized in the extreme case where the 
toroidal magnetic field B, attains a maximum possible value ~ Bp due to magnetic turbulent dif- 
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fusion. Usually, the coupling of matter with the magnetosphere is mediated by different plasma 
instabilities whose characteristic times are too short for substantial toroidal field growth. As we 
discussed above in Section [ITl the shell is very hot near the bottom, so without cooling at the mag- 
netospheric boundary it is marginally stable with respect to the interchange instability, according 
to the calculations by Arons and Lea (1976a). Due to Compton cooling by X-ray emission from the 
neutron star poles, plasma enters the magnetosphere with a characteristic time-scale determined 
by the instability increment. Since the most likely instability is the Rayleigh-Taylor instability, 
the growth rate scales as the Keplerian angular frequency ojk = ^[GMJW. The time- scale of the 
instability can be normalized as 

hnst — Tn~\ ' (22) 

The toroidal magnetic field increases with time as 

B, = KmBp{oj^-oj*)tin,,, (23) 

where (6) is a numerical coefficient which takes into account the degree and angular dependence 
of the coupling of matter with the magnetosphere in which the angle between the neutron star spin 
axis and the magnetic dipole axis is 6. Then, the neutron star spin frequency change in this regime 
reads 

IcU* = KmK2^^^^^^ . (24) 
Using the definition of Ra [Eq. (fTTI) ! and cok, the spin-down formula can be recast to the form 
I(b* = ZMR\{oj,r, - of) (25) 
Here the dimensionless coefficient Z is 

Z = — — (/r(y, mt){\ + ym, ) . (26) 

f{u) Ay 

Taking into account that the matter falling onto the neutron star surface brings angular mo- 
mentum zMR\of (see Eq. (fTSi) above), we arrive at 

Ico* = ZMRl{o)„, - of) + zMRlof (27) 

Clearly, to remove angular momentum from the neutron star through this kind of a static shell, 
Z must be larger than z- Then the neutron star can spin-down episodically (we shall precise this 
statement below). Oppositely, if Z < z, the neutron star can only spin up. 

When a hot shell cannot be formed (at high accretion rates or small relative wind velocities, 
see e.g. Sunyaev (1978)), free-fall Bondi accretion with low angular momentum is realized. No 
angular momentum can be removed from the neutron star magnetosphere. Then Z = z and Eq. (1271) 
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takes the simple form I6f = ZMR\ 



(jj,„, and the neutron star in this regime can spin-up up to 



about (jjKiRA) independent of the sign of the difference of the angular frequencies co,,, - of at the 
magnetopsheric boundary. Due to conservation of specific angular momentum, o),,, = tOBiRB/RA)^, 
so in this case the spin evolution of the NS is described by equation 



where Z plays the role of the specific angular momentum of the captured matter. For example, 
in early work by lUarionov & Sunyaev 1975 Z ^ 1/4. However, detailed numerical simulations 
of Bondi-Littleton accretion in 2D (e.g. Fryxell and Taam, 1988, Ho et al. 1989) and in 3D (e.g. 
RufFert, 1997, 1999) revealed that due to inhomogeneities in the incoming flow, a non-stationary 
regime with an alternating sign of the captured matter angular momentum can be realized. So the 
sign of Z can be negative as well, and alternating spin-up/spin-down regimes can be observed. 
Such a scenario is frequently invoked to explain the observed torque reversals in X-ray pulsars 
(see the discussion in Nelson et al. 1997). We repeat that this could indeed be the case for large 
X-ray luminosities > 4 x 10^^ erg/s when a convective quasi-hydrostatic shell cannot exist due to 
strong Compton cooling near the magnetospheric boundary. 

When a hot shell is formed (at moderate X-ray luminositities below ~ 4 x 10^^ erg/s, see 
Eq. (l59l) below), the angular momentum from the neutron star magnetosphere can be transferred 
away through the shell by turbulent viscosity and convective motions. So we substitute ojm from 
Eq. (O into Eq. ^ to obtain 



This is the main formula which we shall use below. To proceed further, however, we need to 
determine the dimensionless coefficients of this equation. In the next section we shall find the im- 
portant factor f(u) that enters the formulae for both Z and Ra, so the only unknown dimensionless 
parameter of the problem will be the coefficient Ki{6). 

3 APPROXIMATE STRUCTURE OF THE INTERCHANGE INSTABILITY REGION 

The plasma enters the magnetosphere of the slowly rotating neutron star due to the interchange 
instability. The boundary between the plasma and the magnetosphere is stable at high temperatures 
T > Tcr, but becomes unstable at T < r„-, and remains in a neutral equilibrium at T = T^- (Eisner 
and Lamb, 1976). The critical temperature is: 



Ico* = ZMcOfiRl 



(28) 




(29) 
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1 COSY Ll,„GM 

nT,r = ^ -— (30) 

2(1 +rm?) kRa Ra 

Here k is the local curvature of the magnetosphere, x is the angle the outer normal makes with 
the radius-vector at a given point, and the contribution of turbulent pulsations in the plasma to the 
total pressure is taken into account by factor (1 + ymj). The effective gravity acceleration can be 
written as 

GM I T 



(31) 

The temperature in the quasi-static shell is given by Eq. (Hj), so the condition for the magnetosphere 
instability can then be rewritten as: 

Tf- = Hr^mt) < 1 . (32) 

According to Arons and Lea (1976a), when the external gas pressure decreases with radius as 
P ~ R~^^^, the form of the magnetosphere far from the polar cusp can be described to within 10% 
accuracy as (cos ^693 (j^gj-g ^ jj^g polar angle counting from the magneto spheric equator). The 
instability first appears near the equator, where the curvature is minimal. Near the equatorial plane 
{A = 0), for a poloidal dependence of the magnetosphere ~ (cos /l)"^^ we get for the curvature 
kpRA = 1-1-1.27. The toroidal field curvature at the magnetospheric equator is fc, = 1. The tangent 
sphere at the equator cannot have a radius larger than the inverse poloidal curvature, therefrom 
kRa = 1.27 at /i = 0. This is somewhat larger than the value of kRa = 7/(2(7 - 1)) = 5/4 = 1.25 ( 
for 7 = 5/3 in the absence of turbulence or for fully isotropic turbulence), but within the accuracy 
limilQ. The contribution from anisotropic turbulence decreases the critical temperature; for exam- 
ple, for 7 = 5/3, in the case of strongly anisotropic turbulence my = 1, = 0, at /I = we obtain 
T/Tcr ~ 2, i.e. anisotropic turbulence increases the stability of the magnetosphere. So initially the 
plasma-magnetospheric boundary is stable, and after cooling to T < the plasma instability sets 
in, starting in the equatorial zone, where the curvature of the magnetospheric surface is minimal. 

Let us consider the development of the interchange instability when cooling (predominantly 
the Compton cooling) is present. The temperature changes as (Kompaneets, 1956, Weymann, 
1965) 

dT T -T, 



dt tc 

where the Compton cooling time is 



(33) 



In Arons and Lea (1976b), the curvature is calculated to be kRa ~ 1.34, still within the accuracy limit 



Quasi-spherical accretion 15 

3 nRlnieC^ , . , 

tc = ^ « \OmRlM-l . (34) 

Here is the electron mass, aj is the Thomson cross section, = O.lMc^ is the X-ray luminos- 
ity, T is the electron temperature (which is equal to ion temperature; the timescale of electron-ion 
energy exchange is the shortest one), is the X-ray temperature and //^ = 0.6 is the molecular 
weight. The photon temperature is = {l/4-)Tcut for a bremsstrahlung spectrum with an expo- 
nential cut-off at Tcut, typically = 3 - 5 keV. The solution of equation Eq. reads: 

T = T, + {T,,--T,)e-'"^ . (35) 

It is seen that for t ~ 2tc the temperature decreases to T^. In the linear approximation and noticing 
that Tcr ~ 30keV :s> ~ 3 keV, the effective gravity acceleration increases linearly with time: 

^^//*-^^- (36) 
Correspondingly, the rate of instability increases with time as 

= J Seffdt = ^-^t ■ (37) 
Let us introduce the mean rate of the instability growth 

<Ui>= ^— = = -—^\- — - . (38) 



t 6 Rl tc 6Rltc \< Ui >, 
Here ^ < 1 and ^Ra is the characteristic scale of the instability that grows with the rate < m, >. So 
for the mean rate of the instability growth in the linear stage we find 

Here we have introduced the free-fall time as 

^3/2 

tff = — ^ . (40) 
^^2GM 

Clearly, later in the non-linear stage the rate of instability growth approaches the free-fall 
velocity. We consider the linear stage first of all, since at this stage the temperature is not too low 
(although the entropy starts decreasing with radius), and it is in this zone that the effective angular 
momentum transfer from the magnetosphere to the shell occurs. At later stages of the instability 
development, the entropy drop is too strong for convection to begin. 

Let us estimate the accuracy of our approximation by retaining the second-order terms in the 
exponent expansion. Then the mean instability growth rate is 



< Ui >- 



etc / 



(41) 



Clearly, the smaller accretion rate, the smaller the ratio tff/tc, and the better our approximation. 
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Now we are in a position to specify the important dimensionless factor f{u): 

UffiRA) 

Substituting Eq. (|39l) and Eq. (I42l) into Eq. (fTTl) . we find for the Alfven radius in this regime: 

. 0.9 X 10^[cm] { -ifl)''" . (43) 

\(r- 1)(1 +rmf)i/r(7,m,) Mie/ 

We stress the difference of the obtained expression for the Alfven radius with the standard one, 
Ra ~ i/^'' I M~^l'^ , which is obtained by equating the dynamical pressure of falling gas to the mag- 
netic field pressure; this difference comes from the dependence of f{u) on the magnetic moment 
and mass accretion rate in the settling accretion regime. 

Plugging Eq. (|43] ) into Eq. (|42l) . we obtain an explicit expression for f{u): 

f{u)^Q33y — 1 . (44) 



4 SPIN-UP/SPIN-DOWN TRANSITIONS 

Now let us consider how the spin-up/spin-down torques vary with changing M. We stress again that 
we consider the shell through which matter accretes to be essentially subsonic. It is the leakage of 
matter through the magneto spheric boundary that in fact determines the actual accretion rate onto 
the neutron star. This is mostly dependent on the density at the bottom of the shell. On the other 
hand, the density structure of the shell is directly related to the density of captured matter at the 
bow shock region, so density variations downstream the shock are rapidly translated into density 
variations near the magnetopsheric boundary. This means that the actual accretion rate variations 
must be essentially independent (for circular or low-eccentricity orbits) of the orbital phase but 
are mostly dependent on variations of the wind density. In contrast, possible changes in Rb (for 
example, due to wind velocity variations or to the orbital motion of the neutron star) do not affect 
the accretion rate through the shell but strongly affect the value of the spin-up torque (see Eq. (|29l)). 
Eq. (|29l ) can be rewritten in the form 

lio* = AM^ - BM^^^^ . (45) 

For the fiducial value Mie = M/IO^^ g/s, the accretion-rate independent coefficients are (in COS 
units) 

13-611 

^ u~v~^"i — ] ( 

(l+(5/3)m?)<A(5/3,m,); "^^^ ' llOd/ 



A « 5325xl0^\0.034)'^-"Ki{e)(b6"{l+i5/3)m^) 
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5 = 5.4 X \0'\\ - zizwm I X < ffSr) (47) 

\(1 + (5/3)m2)i/r(5/3,m,)/ '^'\100s/ 

(from now on in numerical estimates we assume y = 5/3). The dimensionless factor 6 < I takes 
into account the actual location of the gravitaional capture radius, which is smaller than the Bondi 
value for a cold stellar wind (Hunt, 1971). This radius can also be smaller due to radiative heating 
of the stellar wind by X-ray emission from the neutron star surface (see below). In the numerical 
coefficients we have used the expression for Z with account for Eq. (|44l) . and Eq. (|43]) for the 
Alfvenic radius. 

Below we shall consider the case Z - z > 0, i.e. B > 0, otherwize only spin-up of the neutron 
star is possible. 

First of all, we note that the function to*(M) reaches minimum at some M^-, By differentiating 
Eq. (l45l) with respect to M and equating to zero, we find 



B 



(48) 



A(3 + 2n) 

At M = Mcr the value of d)* reaches an absolute minimum (see Fig. [I]). 
It is convenient to introduce the dimensionless parameter 

and rewrite Eq. (l45l) in the form 



. 3+22 3+2n 

I(b* = AMc^' y— 

■ y 

where the frequency derivative vanishes aty = jq: 



(50) 



yQ = \—; — • (51) 



3 

The qualitative behaviour of co* as a function of y is shown in Fig. \T\ 
Let us then vary Eq. (|50l) with respect to y: 

I{66j*) = = ^^AM,;' y-— \l - j {5y) . (52) 

We see that depending on whether y > 1 or y < 1 , correlated changes of 56f with X-ray flux should 
have different signs. Indeed, for GX 1-1-4 in Gonzalez-Galan et al (2011) a positive correlation 
of the observed 6P with 6M was found using Fermi data. This means that there is a negative 
correlation between 5(jf and 6M, suggesting y < 1 in this source. 
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Figure 1. A schematic plot of (j* as a function ofy [Eq. )50H . In fact, as y — * 0, ui* approaclies some negative w*, since tlie neutron star enters the 
propeller regime at small accretion rates. 



5 DETERMINATION OF THE NEUTRON STAR MAGNETIC FIELD AND OTHER 
PARAMETERS IN THE SETTLING ACCRETION REGIME 

Most X-ray pulsars rotate close to their equilibrium periods, i.e. the average of = 0. Near the 
equilibrium, in the settling accretion regime from Eq. (|45l) we obtain: 



(ecj) 
^30 



l-z/Z 



0.0986 • (0.034)(2-«)a;(l + (5/3)m^)Y'' (PJlOOsY' /^i6(l + (5/3)m2)i/.(5/3,m,)\^ ( yfdY 



V8 J 



Once the magnetic field of the neutron star is estimated for any specific system, we can cal- 
culate the value of the Alfven radius Ra [Eq. (|43]) 1 and the important numerical coefficient f(u) 
[Eq. (|44l) 1. The coupling constant Ki(9) is evaluated from Eq. (|46l ). in which the left-hand side can 
be independently calculated using Eq. (|52| ) measured aty = yo (where d>* = 0): 

'(f)| 



(54) 



/ \ ■ 3+2ii / -(2ji-1) \ 

The coefficient Z is then determined from Eq. (|26|) . The dimensionless factor relating the toroidal 
and poloidal magnetic field is also important. Near the equilibrium we have to^ - co* = -{zlZ)(jf, 
so Eq. (|23]) can be written as 

^ = Kmi-]i^^] ( ] (55) 

Bp ' \z)\a)KiRA)l V2(l +(5/3)m2)i/r(5/3,m,) \^^^(^a)/ ' 

Calculated values for all parameters and coefficients discussed above are listed for specific 
wind-fed pulsars in Table 1 below. 
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5.1 Low-luminosity X-ray pulsars with torque reversal 

Let us consider X-ray pulsars with persistent spin-up and spin-down episodes. We shall assume 
that a convective shell is present during both spin-up (with higher luminosity) and spin-down (with 
lower luminosity), provided that at the spin- up stage the mass accretion rate does not exceed 
as derived above. 

Suppose that in such pulsars we can measure the average value of co*\su and 6f\sci as well as the 
average X-ray flux during spin-up and spin-down, respectively. 

Let at some specific value yi < yo the source be observed to spin-down: 



/a;1,, = AM„." y " 



3+2/1 3+2/1 
1 



1-1^ 



InllU 



At some 3^2 > Jo, the neutron star starts to spin-up: 

, 2«/ll\ 



>0, 3^1 < Jo = I I . (56) 



3+2/1 3+2/1 



< , > 3^0 (57) 



3+2/1 

\X\ = X— 



2n/ll 

(58) 



^,2«/ll 



Using the observed spin-down/spin-up ratio co*\sd/co*\.nt = X and the corresponding X-ray lumi- 
nosity ratio Lx{sd) I Lx{su) = yilyi = x < find by dividing Eq. (l56l) with Eq. (l57l) and after 
substituting 3^2 = y\lx 

Solving this equation, we obtain yi and yj through the observed quantities \X\ and x. We stress that 
so far we have not used the absolute values of the X-ray luminosity (the mass accretion rate), only 
the ratio of the X-ray fluxes during spin-up and spin-down. 

Here we should emphasize an important point. When the accretion rate through the shell ex- 
ceeds some critical value M > M*, the flow near the Alfven surface may become supersonic, a 
free-fall gap appears above the magnetosphere, and the angular momentum can not be removed 
from the magnetosphere. In that case the settling accretion regime is no longer realized, a shock is 
formed in the flow near the magnetosphere (the case studied by Burnard et al. 1982). Depending 
on the character of the inhomogeneities in the captured stellar wind, the specific angular momen- 
tum of the accreting matter can be prograde or retrograde, so alternating spin-up and spin-down 
episodes are possible. Thus the transition from the settling accretion regime (at low X-ray lumi- 
nosities) to Bondi-Littleton accretion (at higher X-ray luminosities) can actually occur before 3' 
reaches yo. Indeed, by assuming a maximum possible value of the dimensioless velocity of mat- 
ter settling /(m)=0.5 (for angular momentum removal from the magnetosphere to be possible, see 
Appendix D for more detail), we find from Eq. (I44l) : 
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M* 3.7 1 I ix''^ . (59) 

\(l+(5/3)m?)^(5/3,m,)/ "^^o 

A similar estimate for the critical mass accretion rate for settling accretion can be obtained 

from a comparison of the characteristic Compton cooling time with the convective time at the 

Alfven radius. 



6 SPECIFIC X-RAY PULSARS 

In this Section, as an illustration of the possible applicability of our model to real sources, we 
consider three particular slowly rotating moderatly luminous X-ray pulsars: GX 301-2, Vela X-1, 
and GX 1-1-4. The first two pulsars are close to the equilibrium rotation of the neutron star, showing 
spin-up/spin-down excursions near the equilibrium frequency (apart from the spin-up/spin-down 
jumps, which may be, we think, due to episodic switch-ons of the strong coupling regime when 
the toroidal magnetic field component becomes comparable to the poloidal one, see Section [2.4.1l) . 
The third one, GX 1-1-4, is a typical example of a pulsar displaying long-term spin-up/spin-down 
episodes. During the last 30 years, it shows a steady spin-down with luminosity-(anti)correlated 
frequency fluctuations (see Gonzalez-Galan et al. 2011 for a more detailed discussion). Clearly, 
this pulsar can not be considered to be in equilibrium. 

6.1 GX 301-2 

GX301-2 (also known as 4U 1223-62) is a high-mass X-ray binary, consisting of a neutron star 
and an early type B optical companion with mass - 40Mo and radius =^ 60Rq. The binary period 
is 41.5 days (Koh et al. 1997). The neutron star is a ~ 680 s X-ray pulsar (White et al. 1976), 
accreting from the strong wind of its companion {Mioss ~ 10"^Mo/yr, Kaper et al. 2006). The 
photospheric escape velocity of the wind is Vesc ~ 500 km/s. The semi-major axis of the binary 
system is a ~ IIORq and the orbital eccentricity e ~ 0.46. The wind terminal velocity was found 
by Kaper et al. (2006) to be about 300 km/s, smaller than the photospheric escape velocity. 

GX 301-2 shows strong short-term pulse period variability, which, as in many other wind- 
accreting pulsars, can be well described by a random walk model (deKool & Anzer 1993). Earlier 
observations between 1975 and 1984 showed a period of ~ 700s while in 1984 the source started 
to spin up (Nagase 1989). BATSE observed two rapid spin-up episodes, each lasting about 30 
days, and it was suggested that the long-term spin-up trend may have been due entirely to similar 
brief episodes as virtually no net changes in the frequency were found on long time-scales (Koh et 
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Figure 2. Torque-luminosity correlation in GX 301-2, di' as a function of BATSE data (20-40 keV pulsed flux) near the equilibrium frequency, see 
Doroshenko et al. (2010). The assumed X-ray flux at equilibrium (in terms of the dimensionless parameter v) is also shown by the vertical dotted 
line. 



al. 1997; Bildsten et al. 1997). The almost 10 years of spin-up were followed by a reversal of spin 
in 1993 (Pravdo & Ghosh 2001) after which the source has been continuously spinning down (La 
Barbera et al. 2005; Kreykenbohm et al. 2004, Doroshenko et al. 2010). Rapid spin-up episodes 
sometimes appear in the Fermi/GBM data on top of the long-term spin-down trend (Finger et al. 
201 1). It can not be excluded that these rapid spin-up episodes, as well as the ones observed in the 
BATSE data, reflect a temporary entrance into the strong coupling regime, as discussed in Section 
2.4.1. Cyclotron line measurements (La Barbera et al. 2005) yield the magnetic field estimate near 
the neutron star surface Bq « 4.4 x 10'^ G (p = \/2BoRl = 2.2 x 10^° G cm^ for the assumed 
neutron star radius Rq = 10 km). 

In Fig. [2] we have plotted to* as a function of the observed pulsed flux (20-40 keV) according 
to BATSE data (see Doroshenko et al. 2010 for more detail). To obtain the magnetic field estimate 
and other parameters (first of all, the coefficient A, see Eq. (|54l )) as described above in Section [51 
we need to know the value of M and the derivative dof IdM or dof Idy. The estimate of M can be 
inferred from the X-ray flux provided the distance to the source is known, and generally this is a 
major uncertainty. We shall assume that near equilibrium a hot quasi-spherical shell exists in this 
pulsar, i.e. the accretion rate is 3 x 10^^ g/s, i.e. not higher than the critical value ^ 4 x 10'^ g/s 
[Eq. (l59l) 1 . While the absolute value of the mass accretion rate is necessary to estimate the magnetic 
field according to Eq. (|53l) (however, the dependence is rather weak, ~ M^^^), the derivative dof Idy 
can be derived from the co* - X-ray flux plot, since in the first approximation the accretion rate 



22 A^. Shakura et al. 



1 

— I — 

+^ 



^— 

I I I I I I L 

0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 

BATSE flux, cts/s 

Figure 3. The same as in Fig.|2]for Vela X-1 (Doroshenko 201 1, private communication). 

is proportional to the observed pulsed X-ray flux. Near the equilibrium (the torque reversal point 
with to* = 0), we find from a linear fit in Fig.^doj* / dy « 4 x 10"^^ rad/s^. 

The obtained parameters (pi, Z, Ki{theta), etc.) for this pulsar are listed in Table 1. We note 
that the magnetic field estimate resulting from our model for n = 2 (boldfaced in Table 1) is fairly 
close to the value inferred from the cylcotron line measurements. We also note that for the case 
n = 3/2 the coupling constants Z and Ki{0) turn out to be unrealistically large and the derived 
magentic field is very small, suggesting that assuming anisotropic turbulence is more realistic than 
using isotropic turbulence with the viscosity described by the Prandtl law. 

6.2 Vela X-1 

Vela X-1 (=4U 0900-40) is the brightest persistent accretion-powered pulsar in the 20-50 keV 
energy band with an average luminosity of Lv « 4 x lO^^erg/s (Nagase 1989). It consists of a 
massive neutron star (1.88 Mq, Quaintrell et al. 2003) and the BO.SIb super giant HD 77581, 
which eclipses the neutron star every orbital cycle of ~ 8.964 d (van Kerkwijk et al. 1995). The 
neutron star was discovered as an X-ray pulsar with a spin period of ~283 s (Rappaport 1975), 
which has remained almost constant since the discovery of the source. The optical companion 
has a mass and radius of ~ 23 Mq and ~ 30 respectively (van Kerkwijk et al. 1995). The 
photospheric escape velocity is Vesc * 540 km/s. The orbital separation is a « 50Rq and the orbital 
eccentricity e « 0.1. The primary almost fills its Roche lobe (as also evidenced by the presence 
of elliptical variations in the optical light curve, Bochkarev et al. (1975)). The mass-loss rate from 
the primary star is 10"^ MJyx (Nagase et al. 1986) via a fast wind with a terminal velocity of ~ 
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1 100 km/s (Watanabe et al. 2006), which is typical for this class. Despite the fact that the terminal 
velocity of the wind is rather large, the compactness of the system makes it impossible for the 
wind to reach this velocity before interacting with the neutron star, so the relative velocity of the 
wind with respect to neutron star is rather low, ~ 700 km/s. 

Cyclotron line measurements (Staubert 2003) yields the magnetic field estimate Bq x 3x10^^ G 
(// = 1.5 X 10^" G cm^ for the assumed neutron star radius 10 km). We shall assume that in this 
pulsar M ^ 3x 10^^ g/s (again for the existence of the shell to be possible). In Fig.[3]we have plotted 
d>* as a function of the observed pulsed flux (20-40 keV) according to BATSE data (Doroshenko 
2011, private communication). As in the case of GX 301-2, from a linear fit we find at the spin- 
up/spin-down transition point X 5.5 X 10"^^ rad/s^. 

The obtained parameters for Vela X-1 are listed in Table 1. As in the case of GX 1-1-4, the 
magnetic field estimate given by our model for an almost iso-angular-momentum rotation law 
(n = 2, boldfaced in Table 1) is close to the value inferred from the cyclotron line measurements. 



6.3 GX 1+4 

GX 1-1-4 was the first source to be identified as a symbiotic binary containing a neutron star (David- 
sen, Malina & Bowyer 1977). The pulse period is ~ 140 s and the donor is an Mill giant (Davidsen 
et al. 1977). The system has an orbital period of 1 161 days (Hinkle et al. 2006) making it the widest 
known LMXB by at least one order of magnitude. The donor is far from filling its Roche lobe and 
accretion onto the neutron star is by capture of the stellar wind of the companion. 

The system has a very interesting spin history. During the 1970's it was spinning up at the 
fastest rate (cosu ~ 3.8 • 10"^' rad/s) among the known X-ray pulsars at the time (e.g. Nagase 1989). 
After several years of non-detections in the early 1980's, it reappeared again, now spinning down at 
a rate similar in magnitude to that of the previous spin-up. This spin-reversal has been interpreted 
in terms of a retrograde accretion disc forming in the system (Makishima et al. 1988, Dotani et 
al. 1989, Chakrabarty et al. 1997). A detailed spin-down history of the source is discussed in the 
recent paper by Gonzalez- Galan et al. (2011). Using our model this behavior can, however, be 
readily explained by quasi-spherical accretion. 

As the pulsar in GX 1-1-4 is not in equilibrium, we cannot directly use our method to estimate 
the magnetic field of the neutron star as described in Section |5l However, as GX 1-1-4 is currently 
experiencing a long-term spin-down trend, the first (spin-up) term in Eq. (|45l) must be smaller than 
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Figure 4. Pulsar frequency (upper panel), deviation of the frequency from the linear fit (middle panel) and pulsed flux in GX 1 +4 from Fermi GBM 
data (M. Finger, piivate communication; see also Gonzalez-Galan et al (2011)) 



the second (spin-down) term, which yields a lower limit on the value of the magnetic field (see 
Table 1). 

Let us use Eq. (l52l) to quantitatively explain the observed anti-c orrelation between the pulsa r 



frequency fluctuations and the X-ray flux, as observed in GX 1-1-4 (IGonzalez-Galan et al, 
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We use a fragment of four-day average Fermi/GBM data on GX 1-1-4 (M. Finger, private com- 
munication; see also Gonzalez-Galan et al (201 1)). The pulsar is currently observed at the steady 
spin-down stage with a mean iOsd « -2.34 x 10 " rad/s (the upper panel). In the middle panel 
of Fig. m deviations from the linear fit are shown. Note that the frequency excursions around the 
mean value is a few microseconds, while the pulsar frequency is much higher, a few milliseconds. 
This means that the frequency derivative is negative at all points within the time interval shown, 
i.e. no occasional spin-ups were observed, even at the highest X-ray flux levels. 

Specifically, let us consider the prominent pulsar frequency change observed between MJD 
55100 and MJD 55200 for around At = 80 days. During this time period, the frequency of 
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the pulsar decreased by A(y* « -3.6 x 10"^ rad (see Fig. H]). From here we find 56f{obs) = 
ls.(jf I ls.t ~ -5.2 X 10~'^ rad/s. Thus, the observed fractional change in the pulsar spin-down rate is 
{6co*IW\Us « -0.2. 

On the other hand, by dividing Eq. (|52l ) with Eq. (|50l) we find the expected relative fluctuations 
of of at a given mean accretion rate (or dimensionless X-ray luminosity y): 



66j* 3 + 2n6M 



W\ theor 11 M 



. ln-\ 



i-(?r 



(60) 



From Fig. |4] (bottom panel) the range of the fluctuations relative to the mean value {SMjM) ^ 
(0.6 - 0.2)/0.4 ^ 1. Then, for the assumed n = 2, the expected amplitude of the relative frequency 
derivative fluctuation would match the observed value if y ^ 0.2 - 0.3. Important is that we find 
y < 1, as must be the case for the observed negative sign of the torque-luminosity correlations. 

Note here that the drop in flux-level during about 20 days at around MJD 55140-55160 should 
be translated to a decrease in to*, which is clearly seen in the upper panel of Fig.lH A more detailed 
analysis using the entire Fermi data-set should be performed. 

Further, note that the short-term spin-up episodes, sometimes observed on top of the steady 
spin-down behaviour (at about MJD 49700, see Fig. 2 in Chakrabarty et al. (1997)) are correlated 
with an enhancement of the X-ray flux, in contrast to the negative frequency-flux correlation dis- 
cussed above. During these short spin-ups, d)* is about half the average 6c»*„ observed during the 
steady spin-up state of GX 1-1-4. The X-ray luminosity during these episodic spin-ups is approxi- 
mately five times larger than the mean X-ray luminosity during the steady spin-down. We remind 
the reader that once M > M,, a free-fall gap appears above the magnetosphere, and the neutron star 
can only spin up. When the X-ray flux drops again, the settling accretion regime is reestablished 
and the neutron star resumes its spinning-down. 



7 DISCUSSION 

7.1 Physical conditions inside the shell 

For an accretion shell to be formed around the neutron star magnetosphere it is necessary that the 
matter crossing the bow shock does not cool down too rapidly and thus starts to fall freely. This 
means that the radiation cooling time tcooi must be longer than the characteristic time of plasma 
motion. 

The plasma heats up in the strong shock to the temperature 
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Table 1. Parameters for the pulsars in Section 6. References for the observed spin periods, binary periods and spin down rate (for GX 1 +4) are 
given in the text as well as discussions on the values used for the wind velocities with respect to the neutron star. The parameters Z, Ki(0) and 
f(u) are defined in Section 2.4.2. Numerical estimates are given for dimensionless parameters 5 = 1, f = 1, tj = 1, y = 5/3 and without turbulence 
(m, = 0). The numbers in boldface are the preferred values for a near iso-angular-momentum rotation law. 



Pulsar GX301 - 2 VelaX - 1 GXl +4 

Measured parameters 



P.(s) 680 283 140 

PB(d) 41.5 8.96 1161 

iv(km/s) 300 700 200 



fly l>0 


4- 


io-'3 


5.5 


10-13 




Mi6(M/10'«) 




3 




3 


1 


Derived parameters 




n = 2 


n = 3/2 


n = 2 


n = 3/2 


n = 2 « = 3/2 




2.7 


0.1 


1.8 


0.16 


> 1.17 > 0.02 


/(«) 


0.42 


0.57 


0.43 


0.54 




Kim 


39 


3700 


36 


1150 




Z 


13 


910 


12 


300 




B,/B„ 


0.1 


0.01 


0.2 


0.03 




RA(cm) 


2 10' 


3- 10** 


1.6 10' 


4.2- 10** 






0.06 


0.004 


0.1 


0.01 





The radiative cooling time of the plasma is 
3kT 



^cool 



(61) 



(62) 



where p is the plasma density, = Y^p/nip is the electron number density ( = 0.6 and « 0.8 
for fully ionized plasma with solar abundance); A is the cooling function which can be approxi- 
mated as 

o,r < lO^K 

i.ox io-24ro55 lo^K < r < lo^K 

6.2 X 10-197-0 6^ 10^ K < r < 4 X 10^ K 
2.5xl0-27r0-5,r>4xl0^K 
(Raymond, Cox & Smith 1976; Cowie, McKee & Ostriker 1981). 

Compton cooling becomes effective from the radius where the gas temperature T, determined 
by the hydrostatic formula Eq. dU, is lower than the X-ray Compton temperature T^. The Compton 
cooling time (see Eq. (I34l) ) is: 



(63) 



tc ~ 1060[s]M 



16 



R 



10.«cm/ ■ 

Above the radius where = T, Compton heating dominates. Taking the actual temperature 
close to the adiabatic one [Eq. dH)], we find ^ 2x 10'° cm. We note that both the Compton and 
photoionization heating processes are controlled by the photoionization parameter ^ (Tarter et al. 
1969, Hatchett et al. 1976) 
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^=-^- (65) 
In most part of the accretion flux, n ~ R~^l^, so ^ ~ Rr^l'^ and independent of the X-ray luminosity 
through the mass continuity equation. For characteristic values we find: 

^«5xlOV(")<^'- (66) 

If Compton processes were eff"ective everywhere, this high value of the parameter ^ would imply 
that the plasma is Compton-heated up to keV-temperatures out to very large distances ~ 10'^ cm. 
However, at large distances the Compton heating time becomes longer than the characteristic time 
of gas accretion: 

;^ = ^£^-20/MMr.'^:«, (67) 

^accr ^ 

which shows that Compton heating is ineff"ective. The gas temperature is determined by photoion- 
ization heating only and the gas can only be heated up to T„y^,x « 5 x 10^ K (Tarter et al. 1969), 
which is substantially lower than ~ 3 keV. The sound velocity corresponding to T^ax is approx- 
imately 80 km/s. 

The eff"ective gravitational capture radius corresponding to the sound velocity of the gas in the 
photoionization-heated zone is 

Rn = — ^ = ~ 3.5 X 10 cm T— . (68) 

Everywhere up to the bow shock photoionization keeps the temperature at a value =^ T^ax- 

If the stellar wind velocity exceeds 80 km/s, a standard bow shock is formed at the Bondi 
radius with a post-shock temperature given by Eq. (|6T1) . If the stellar wind velocity is lower than 
this value, the shock disappears and quasi-spherical accretion occurs fromi?^. 
The photoionization heating time at the eff"ective Bondi radius 3 x 10^^ cm is 

tpi ~ ^- ~ 2 X 10 [s]Mi6 . (69) 

(here hveff ~ 10 keV is the characteristic photon energy, ^ is the eff'ective photoionization po- 
tential, aeff ~ 10"^"* cm^ is the typical photoionization cross-section, Uy = L/{4nR^hVeffC) is the 
photon number density). The photoionization to accretion time ratio at the effective Bondi radius 
is then 

« Omf(u)M^^ . (70) 

''accr 

At wind velocities > 80 km/s the bow shock stands at the classical Bondi radius Rb inside 
the eff'ective Bondi radius R*^ determined by Eq. (|68l) . The cooling time of the shocked plasma at 
Rb expressed through the wind velocity v„ is: 
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t.ooi « 4.7 X . (71) 

The photoionization heating time in the post-shock region can also be expressed through the stellar 
wind velocity: 

« 3.5 X 10^[s]M->7^ (72) 

The comparison of these two timescales implies that at low velocities radiative cooling is important 
and the regime of free-fall accretion with conservation of specific angular momentum is realized. 

So, at low wind velocities the plasma cools down and starts to fall freely. As the cold plasma 
approaches the gravitating center, photoionization heating becomes important and rapidly heats 
up the plasma to T^ax « 5 x 10^ K. Should this occur at a radius where Tmax < GM/CRR), the 
plasma continues its free fall down to the magnetosphere, still with the temperature T^ax, with 
the subsequent formation of a shock above the magnetosphere. However, if T,„ax is above the 
adiabatic temperatures at this radius, the settling accretion regime will be established even for low 
wind velocities. 

For high-wind stellar velocities v„ > 100 km/s, the post-shock temperature is higher than T^ax, 
photoionization is unimportant, and the settling accretion regime is established if the radiation 
cooling time is longer than the accretion time. From a comparison of these timescales, we find 
the critical accretion rate as a function of of the wind velocity below which the settling accretion 
regime is possible: 

M**<0.12v^^ (73) 

Here we stress the difference of the critical acccretion rate M** from M* derived earlier. At 
M > M**, the plasma rapidly cools down in the gravitational capture region and free-fall accretion 
begins (unless photoionization heats up the plasma above the adiabatic value at some radius), while 
at M > M* ^ 4 X 10'^ g/s determined by Eq. (l59l) a free-fall gap appears immediately above the 
neutron star magnetosphere. 



7.2 On the possibility of the propeller regime 

The very slow rotation of the neutron stars in GX 1+4, GX 301-2 and Vela X-1 (a)*(RA) < a;^(i?A)) 
makes it hard to establish the propeller regime where matter is ejected with parabolic velocities 
from the magnetosphere during spin-down episodes. 

Let us therefore start with estimating the important ratio of viscous tensions (~ BtBp) to the gas 
pressure (~ B^) at the magneto spheric boundary. This ratio is proportional to B,IBp (see Eq. (1551) ) 
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and is always much smaller than 1 (see Table 1), i.e. only large-scale convective motions with the 
characteristic hierarchy of eddies scaled with radius can be established in the shell. When of > 
(jJk{Ra), the propeller regime (without accretion) must set in. In that case the maximum possible 

braking torque is fJ-^/R\ strong coupling between the plasma and the magnetic field. 

Note that in the propeller state, interaction of the plasma with the magnetic field is in the strong 
coupling regime, i.e. where the toroidal magnetic field component Bt is comparable to the poloidal 
one Bp. It can not be excluded that a hot iso-angular-momentum envelope could exist in this case 
as well, which would then remove angular momentum from the rotating magneto sphere. If the 
characteristic cooling time of the gas in the envelope is short in comparison to the falling time of 
matter, the shell disappears and one can expect the formation of a 'storaging' thin Keplerian disc 
around the neutron star magnetosphere (Sunyaev & Shakura 1977). There is no accretion of matter 
through such a disc. It only serves to remove angular momentum from the magnetosphere. 

7.3 Effects of the hot shell on the X-ray energy and power spectrum 

The spectra of X-ray pulsars are dominated by emission generated in the accretion column. The 
hot optically thin shell produces its own thermal emission, but even if all gravitational energy were 
released in the shell, the ratio of the X-ray luminosity from the shell to that of the accretion column 
would be about the ratio of the magnetosphere radius to the NS radius, i.e. one percent or less. In 
reality, it is much smaller. The shell should scatter X-ray radiation from the accretion column, but 
for this effect to be substantial, the Comptonization parameter y must be of the order of one. The 
Thomson depth in the shell is, however, very small. Indeed, from the mass continuity equation and 
Eq. (|43]) for the Alfven radius and Eq. (|44l) for the factor f{u), we get: 

Tt = n,(R)crTdR « 3.2 x lO'^Mj^^ V^^^^ • 
Jra 

Therefore, for the temperature near the magnetosphere [Eq. dH)] the parameter y is 

nieC^ 

This means that the X-ray spectrum of the accretion column should not be significantly affected 
by scattering in the hot shell. 

The large-scale convective motions in the shell introduce an intrinsic time-scale of the order of 
the free-fall time that could give rise to features (e.g. QPOs) in the power spectrum of variability. 
QPOs were reported in some X-ray pulsars (see Marykutty et al. 2010 and references therein). 
However, the expected frequency of the QPOs arising in our model would be of the order of mHz, 
much shorter than those reported. 
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A stronger effect can be the appearance of a dynamical instability of the shell on this time scale 
due to increased Compton cooling and hence increased mass accretion rate in the shell. This may 
result in a complete collapse of the shell resulting in an X-ray outburst with duration similar to the 
free-fall time scale of the shell (~ 1000 s). Such a transient behaviour is observed in supergiant fast 
X-ray transients (SFXTs) (see Ducci et al. 2010). This interesting issue depends on the specific 
parameters of the shell and needs to be further investigated. 

7.4 Can accretion discs (prograde or retrograde) be present in these pulsars? 

The analysis of real pulsars carried out earlier suggested that in a convective shell an iso-angular- 
momentum distribution is the most plausible. Therefore, we shall below consider only this case, i.e. 
the rotation law co ~ R~^. As follows from Eq. (|29l) . at oj* = the equilibrium angular frequency 
of the neutron star is 

1 IRb^^ 



We stress that such an equilibrium in our model is possible only when a shell is present. At high 
accretion rates M > M* - 4 x 10'^ g/s accretion proceeds in the free-fall regime (with no shell 
present). 

Using Eq. (|53l) . the equlibrium period for quasi-spherical settling accretion can be recasted to 
the form 

10d/\(l+(5/3)m2),A(5/3,m,)MiJ \V^j a>(l + (5/3)m2) ' 
For standard disc accretion, the equilibrium period is 



P. « 1000[s],jr' ...... ... 1 (^1 ^^-l'^ . • (75) 



P,,,,^lsfil'^'M;^" , (76) 

and the long periods observed in some X-ray pulsars can be explained assuming a very high 
magnetic field of the neutron star. Retrograde accretion discs are also discussed in the litera- 
ture (see, e.g.. Nelson et al. (1997) and references therein). Torque reversals produced by pro- 
grade/retrograde discs can in principle lead to very long periods for X-ray pulsars even with stan- 
dard magnetic fields. Retrograde discs can be formed due to inhomogeneities in the captured stellar 
wind (Ruff"ert 1997, 1999). This might be the case at high accretion rates when hot quasi- spherical 
shell cannot exist. In the case of GX 1-1-4, however, it is highly unlikely to observe a retrograde 
disk on a time scale much longer than the orbital period (see a more detailed discussion of this 
issue in Gozalez-Galan et al. (201 1)). In the case of GX 301-2 and Vela X-1, the observed positive 
torque-luminosity correlation (see Figs. [2] and [3]) rules out a retrograde disc as well. 
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To conclude the discussion, we should mention that real systems (including those considered 
here) demonstrate a complex quasi-stationary behaviour with dips, outbursts, etc. These considera- 
tions are beyond the scope of this paper and definitely deserve further observational and theoretical 
studies. 

8 CONCLUSIONS 

In this paper we have presented a theoretical model for quasi-spherical subsonic accretion onto 
slowly rotating magnetized neutron stars. In this model the accreting matter is gravitationally cap- 
tured from the stellar wind of the optical companion and subsonically settles down onto the rotating 
magnetosphere forming an extended quasi-static shell. This shell mediates the angular momen- 
tum removal from the rotating neutron star magnetosphere during spin-down states by large-scale 
convective motions. A detailed analysis and comparison with observations of two specific X-ray 
pulsars GX 301-2 and Vela X-1 demonstrating torque-luminosity correlations near the equilib- 
rium neutron star spin period shows that most likely strongly anisotropic convective motions are 
established, with an almost iso-angular-momentum distribution of rotational velocities a> ~ R~^. A 
statistical analysis of long-period X-ray pulsars with Be-components in SMC (Chashkina & Popov 
2011) also favored the rotation law oj ~ R~^. The accretion rate through the shell is determined 
by the ability of the plasma to enter the magnetosphere. The settling regime of accretion which 
allows angular momentum removal from the neutron star magnetosphere can be realized for mod- 
erate accretion rates M < ^ 4 x 10^^ g/s. At higher accretion rates a free-fall gap above the 
neutron star magnetosphere appears due to rapid Compton cooling, and accretion becomes highly 
non-stationary. 

From observations of the spin-up/spin-down rates (the angular rotation frequency derivative 
to*, or dco* IdM near the torque reversal) of long-period X-ray pulsars with known orbital periods 
it is possible to determine the main dimensionless parameters of the model, as well as to estimate 
the magnetic field of the neutron star. Such an analysis revealed good agreement between the 
magnetic field estimates in the pulsars GX 301-2 and Vela X-1 obtained using our model and 
derived from the cyclotron line measurements. 

In our model, long-term spin-up/spin-down as observed in some X-ray pulsars can be quan- 
titatively explained by a change in the mean mass accretion rate onto the neutron star (and the 
corresponding mean X-ray luminosity). Clearly, these changes are related to the stellar wind prop- 
erties. 
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The model also predicts the specific behaviour of the variations in 56/ , observed on top of 
a steady spin-up or spin-down, as a function of mass accretion rate fluctuations 5M. There is 
a critical accretion rate Mcr below which an anti-correlation of 6af with 6M should occur (the 
case of GX 1+4 at the steady spin-down state currently observed), and above which Saf should 
correlate with 5M fluctuations (the case of Vela X-1, GX 301-2, and GX 1+4 in the steady spin- 
up state). The model explains quantitatively the relative amplitude and the sign of the observed 
frequency fluctuations in GX 1+4. 



APPENDIX A: THE STRUCTURE OF A QUASI-SPEHRICAL ACCRETING 
ROTATING SHELL 

Al Basic equations 

Let us first write down the Navier-Stokes equations in spherical coordinates R, 6, <p. Due to the 
huge Reynolds numbers in the shell (~ 10^^ - 10^^ for a typical accretion rate of 10" g/s and mag- 
netospheric radius ~ 10^ cm), there must be turbulence. In this case the Navier-Stokes equations 
are usually called the Reynolds equations. In the general case, the turbulent viscosity can depend 
on the coordinates, so the equations take the form: 

1. Mass continuity equation: 

dt R^dR^ ^ ' RsinGde ^ RsinO d<p 

2. The i?-component of the momentum equation: 

2 2 

duR duR ueduR duR u^ + Ug gM 

ot oR R 89 Rsmd d(f> R R^ 

3. The ^-component of the momentum equation: 

dug dug UgdUg Us dUg URUg - ulcOtO 

ot oR R 89 R sm 6 ocj) R 

4. The ^-component of the momentum equation: 

dus dus ugdud, Us dus UrUs + u^Ug cot 6 

-r- + Ur-t-^ + — -r-^ + — T-^ + = A^^ (A4) 

dt dR R de RsinO d(f) R ^ 

Here the force components (including viscous force and gas pressure gradients) read: 



— T — — in vvrr] -i — K yy RG sill u) -\ — — Wrs (A5) 

R^dR^ ' sinORde sm9Rd(p ^ R R ^ ^ 



d „. Wgg _ Ws 

.Wg, 
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The components of the stress tensor include a contribution from both the gas pressure Pg 
(assumed to be isotropic) and the turbulent pressure P' (generally anisotropic). In their definition 
we shall follow the classical treatment by Landau and Lifshitz (1986), but with the inclusion of 
anisotropic turbulent pressure: 

Wrr = -Pg - + 2pv,-^ - -py,divu (A8) 

, / 1 dUg Ur\ 2 

Wee = -Pg - P'ee + ^P^t + " ^P^Aiyn (A9) 

117 (IduR due ue\ 

\R smd d(f> R 36 R j 
w / 1 duR du^ uA 

Wr^= pVt\-^-T—— + ^ - — \ (A13) 



^Rsine d4> dR R 

In our problem the anistropy of the turbulence is such that PJ^^ = i^ii, P'^g = P'^^ = P'j^. The 

turbulent pressure components can be expressed through turbulent Mach numbers and will be 

given in Appendix D. 

divu in spherical coordinates is: 

1 (9 / T \ Id 1 dUd, 

R^ oR^ ' R sm BoQ R sm Q 0(p 

A2 Symmetries of the problem 

d d 

We shall consider axially-symmetric ( — = 0), stationary ( — = 0), and only radial accretion 

d(p ot 

{ue = 0). Under these conditions, from the continuity equation Eq. (|A1I) we obtain: 

M = AnR^puR = const . (A15) 

The constant here is determined from the condition of plasma leakage through the magnetosphere. 

Let us rewrite the Reynolds equations under the above assumptions. The i?-component of the 
momentum Eq. (IA2I) equation becomes: 



P 



' duR '^l^ 

" R 



-PST + -5^^ (R Wrr) + -T-r^^ (WRe sm9)-— ^ (A16) 

R'^ R^oR^ ' smORoO R R 



The ^-component of the momentum equation: 

M^COt^ 1 (5 / . V 1 d Waa 

-P^^ = (R'WeR) + (Wee sin0) - cot0^ (A17) 

R R^oR^ ' smd R 86 R 
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The ^-component of the momentum equation: 

The components of the stress tensor with anisotropic turbulence take the form: 

1 duR 

^Re=P^'R-^ (A22) 

Generally, there are two separate cases - with isotropic turbulence and anisotropic turbulence, 
which shall be treated differently. First we consider the simplest case of isotropic turbulence with 
Prandtl law for turbulent viscosity. The more general case of anisotropic turbulence will be dis- 
cussed separately in Appendix B. 



A3 Viscosity prescription for isotropic turbulence 

We consider an axisymmetric flow with a very large Reynolds number. By generalizing the Prandtl 
law for the turbulent velocity obtained for plane parallel flows, the turbulent velocity scales as 
Ut ~ ltR(da>/dR). From the similarity laws of gas-dynamics we assume It ~ R, so 



u, = CiR^ 



dco 



dR 



(A25) 



We note that in our case the turbulent velocity is determined by convection, so <, 0.5uff (see 
Appendix D). This implies that the constant Ci scales as 



Ci ~ u,l{u^), 

and can be very large since {u^) <^ Uf So the turbulent viscosity coefficient reads: 



do) 



dR 



(A26) 



(A27) 



Here C2 « 1 /3 is a numerical factor originating from statistical averaging. Below we shall combine 
Ci and C2 into the new coefficient C = C1C2, which can be much larger than one. Due to interaction 
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of convection cells in the sheared flux the amplitude of turbulent motions is thus much larger than 
the average rotation velocity. 



A4 The angular momentum transport equation and the rotation law inside the shell 

A similar problem (that of a rotating sphere in a viscous fluid) was solved in Landau and Lifshitz 
(1986). It was found that the variables are separated and u^{R,6) = U^{R) s'mO. Note that the 
angular velocity co{R) = Utj,(R)/R is independent of 9. Our problem is different from that of the 
sphere in a viscous fluid in several respects: 1) there is a force of gravity present, 2) the turbulent 
viscosity varies with R and can in principle depend on 9, and 3) there is radial motion of matter 
(accretion). These differences lead, as will be shown below, to the radial dependence U^(R) oc 
R~^^^. (We recall that for a rotating sphere in a viscous fluid oc R~^). 

Let us start with solving Eq. (IA18I ). First, we note that for u^(9) ~ sin 9, according to Eq. (IA23I ), 
Wg^ = 0. Next, making use of the continuity equation Eq. (|D8I) and the definition of angular 
velocity, we rewrite Eq. (|A18I) in the form of angular momentum transfer by viscous forces: 

E±^r2 = ^JL±R^w,, (A28) 
RdR R dR ^ 

Now we can integrate Eq. (IA28I) with respect to R to get 

MojR^ = 4nR^WR^ + D , (A29) 

where D is an integration constant. In the case of isotropic turbulence, the viscous stress component 
(force per unit area) from Eq. (IA24I) is 

WR<i,=pytR^ (A30) 
so we obtain 

MojR^ = AnpvtR^^ + D . (A3 1) 

This equation for angular mometum transport by turbulent viscosity is similar to that for disc 
accretion (Shakura and Sunyaev 1973), but different due to spherical symmetry of the problem 
under consideration. 

The left part of Eq. (IA31I) is simply advection of specific angular momentum averaged over the 
sphere ( 1 /2 ojR^ sin^ 9 sin 9d9 = 1 / 3coR^) by the average motion toward the gravitational center 
(accretion). M is negative as well as ||. The first term on the right describes transport of angular 
momentum outwards by turbulent viscous forces. 

The constant D is determined from the equation 
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D = KmK,^'^^^^ (A32) 



(see Eq. (124)) in the text). We consider accretion onto a magnetized neutron star. When D < 0, the 
advection term in the left part of Eq. (IA31I) dominates over viscous angular momentum transfer 
outwards. Oppositely, when D > 0, the viscous term in the right part of Eq. (IA31I) dominates. 
In the case of M = (no plasma enters the magneto sphere), there is only angular momentum 
transport outwards by viscous forces. 

Now let us rewrite Eq. (|A32I) in the form 

^^^^'^'''^^'^^ 

and use the pressure balance condition 

P{Ra) = P,(i?A)(l + W) = = . (A34) 

8;r 27r7?^ 

Using the mass conitnuity equation in the form 
\M\ = AnR^pfiu) ^^GmJr , 

and the expression for the gas pressure Eq. ([8]), we write the integration constant DI\M\ in the form 
D (r-1) {oj,n-of)Rl , 

-^=K, ^-^ '-ipiy, m,)- ^(1 + ym ) . (A35) 

\M\ 7 2 V2/(m) 

Let us consider the case near the neutron star rotation equilibrium to* = 0. In this case accord- 
ing to Eq. (|27|) 

z 

CO,n - CO* = - — CO* , (A36) 

so using definition Z [Eq. (|26l)1. we obtain: 

= -zR\uf . (A37) 
\M\ ^ 

We emphasize that the value of the constant D is fully determined by the dimensionless specific 
angular momentum of matter at the Alfven radius z- 



A5 Angular rotation law 

Let us use Eq. (IA31I) to find the rotation law ooiK). At large distances R ^ R^ (we remind the 
reader that Ra is the bottom radius of the shell) the constant D is small relative to the other terms, 
so we can set D ^ 0. Thus, to obtain the rotation law we shall neglect this constant in the right part 
of Eq. (|A31I) . Next, we substitute Eq. (|A27I) and make use of the solution for the density (which, 
as we shall show below, remains the same as in the hydrostatic solution) 
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p{R) = P{Ra)[-j) (A38) 
in Eq. (IA31I) to obtain 

\m\ ojR^ = AnpiRA) ( CR' i . (A39) 



Rl \dR_ 
After integrating this equation, we find 

4^1/2 

2oj"' = ±-—+D, (A40) 
where 

K= 3, (A41) 

and Di is some integration constant. In Eq. (IA40I) we use only the positive solution (the minus sign 
with constant Di > would correspond to a solution with the angular velocity growing outwards, 
which is possible if the pulsar rotation is zero). If Di 4^ 0, at large R:^ Ra (in the zone close to the 
bow shock) the solid body rotation law would lead to a; ^ const ^ tos- (However, we remind the 
reader that our discussion is not applicable close to the bow shock region.) At small distances from 
the Alfvenic surface the effect of this constant is small and we shall neglect it in the calculations 
below. Then we find 

4 \M\ (Ra\^'^ 

ojiR) = — — (A42) 

9 Anp{RA)CR\\ R I 

i.e. the quasi-Keplerian law a){R) = oj,n(RA/R)^^^- The constant u),„ in the solution given by 
Eq. (IA42I) is obtained after substituting M from the continuity equation at R = Ra into Eq. (IA42I) : 

_ ~ . 4 , \ur(Ra)\ 

(xjm = o)oj{Ra) = T,<^—p^ — • (A43) 

(Here we have introduced the correction factor > 1 to account for the deviation of the exact 
solution from the Keplerian law near Ra). 

As Ur(Ra) is smaller than the free-fall velocity, the above formula implies that a»,„ < (x)k{Ra)-, 
lower than the Keplerian angular frequency. For self-consistency the coefficient C in the Prandtl 
law is determined, according to Eq. (IA43|) . by the ratio of the radial velocity ur to the rotational 
velocity of matter Uf 
^ AJur{Ra)\ AJur{Ra)\ 

C = -(x> = -oj . (A44) 

9 CO,nRA 9 U^{Ra) 

Note that this ratio is independent of the radius R and is actually constant across the shell. Indeed, 
the radial dependence of the velocity ur follows from the continuity equation with account for the 
density distribution Eq. (IA38I) 
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(A45) 



For the quasi-Keplerian law u^{R) ~ l/R^I^, so the ratio ur/U^ is constant. Finally, the angular 
frequency of the shell rotation near the magneto sphere a),„ is related to the angular frequency of 
the motion of matter near the bow-shock as 



In fact, when approaching R^, the integration constant D (which we neglected at large distances 
R » Ra) should be taken into account. The rotational law will thus be different from the Keplerian 
one. 

We stress the principal difference between this regime of accretion and disc accretion. For disc 
accretion the radial velocity is much smaller than the turbulent velocity, and the tangential velocity 
is almost Keplerian and is much larger than the turbulent velocity. The radial velocity in the quasi- 
spherical case is not determined by the rate of the angular momentum removal. It is determined 
only by the "permeability" of the neutron star magnetosphere for infalling matter. In our case we 
assumed it to be of the order of the convection velocity. The tangential velocity for the obtained 
quasi-Keplerian law is much smaller than the velocity of convective motions in the shell. Note also 
that in the case of disc accretion the turbulence can be parametrized by only one dimensionless 
parameter a ~ uj/u^^ with < or < 1 (Shakura & Sunyaev 1973). The matter in an accretion 
disc differentially rotates with a supersonic (almost Keplerian) velocity, while in our case the shell 
rotates differentially with a substantially subsonic velocity at any radius, and the turbulence in the 
shell is essentially subsonic. Of course, our case with an extended shell is strongly different from 
the regime of freely falling matter with a standing shock above the magnetosphere (Arons & Lea 
1976a). 

A6 The case without accretion 

Now let us consider the case when the plasma can not enter the magnetosphere and no accretion 
onto the neutron star occurs. This case is similar to the subsonic propeller regime considered by 
Davies and Pringle (1981). Eq. (|A31|) then takes the form: 



(We remind the reader that the constant D is determined by the spin-down rate of the neutron star, 
D = Id)* < 0). Solving this equation as above, we find for the rotation law in this case: 



_ . (Rb 

iOm — OJCOb — 



3/2 



(A46) 



= 4npv,R^— + D . 



(A47) 
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oj(R) = oj^l^^j , (A48) 



where 

l7rp(RA)^iRA)Ri 
From Eq. (lA27l) we find 

7 
4 

so for a)„ we obtain: 



= ^ ' '„ ■ (A49) 



Vr(RA) = -Coj,„Rl , (A50) 



_ 2 / /|d 



7 \;rCp(i?A)^ 



1/2 



^m = - ^ „ „J . (A51) 



On the other hand, ojm is now related to the bow-shock region parameters as 

OJ,n = COB\^] . (A52) 



R 



APPENDIX B: THE CASE WITH ANISOTROPIC TURBULENCE 

The Prandl rule for viscosity used above that relates the scale and velocity of turbulent pulsations 
with the average angular velocity is commonly used when the turbulence is generated by the 
shear itself. In our problem, the turbulence is initiated by large-scale convective motions in the 
gravitational field. During convection, a strong anisotropic turbulent motions can appear (the radial 
dispersion of chaotic motions could be much larger than the dispersion in the tangential direction), 
and the Prandl law could be inapplicable. Anisotropic turbulence is much more complicated and 
remains a poorly studied case. 

As a first step, we can adopt the empirical law for Wr^ suggested by Wasiutinski (1946): 

do) 

Wr^ = 2p(-Vf + Vr)C0 + VrpR— (B 1) 

uK 

where the radial and tangential kinematic viscosity coefficients are 

Vr = C||<|Mf||)i? 
Vr = CAK\)R 

respectively. Dimensiionless constants Cy and C± are of the order of one. In the isotropic case, 
Vr = V,, Wr^ ~ doj/dR, and in the strongly anisotropic case, v,- » v,, Wr^ ~ d(a)R^)/dR. Using 
these definitions, let us substitute Eq. (IBll) in Eq. (IA29I ) and rewrite the latter in the form: 

\ M I " M dR \M\ 
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We note that due to self- similarity of the shell structure m|| ~ u'j_ ~ Ur ~ R^^^^, so that the ratios 
(|m[iI)/mr and {\u'^\)/ur are constant. In this case the obvious solution to this equation reads: 



D 



1 



m 1 - 



<KI) 



CObRI + 



D 



1 



\M\ 1 - 2Cx 



<KI) 



|^jC|,(i«Ji> 



l-2Cj_ 



(B3) 



(here the integration constant is defined such that coiRs) = ojb). 

Now let us consider the equilibrium situation where d>* = 0. In this case, as we remember. 



D 2 

--7- = -zofR^ ,iOm = (l- z/Z)a)* . 
\M\ 

First, let us consider the case of strongly anisotropic, almost radial turbulence where (|m'j^|) = 0. 
In this case, the specific angular momentum at the Alfven radius is 



1 + 



C||<|ii,l> 



^ I " _ 1 

Ra 



/J 



C|i<K,l> 



(B4) 



l-z/Z 

It is seen that in the case of very weak accretion (or, in the limit, when the accretion through 
magnetosphere ceases at all), \ur\ «: C||<|m^|), an almost iso-angular-momentum rotation in the 
shell is realized. 

The next case is where the anisotropy is such that Cj_(\u^j^[)/\ur\ = 1/2. Then we have strictly 
iso-angular-momentum distribution: CL>mR\ = (jl>bR\- 

If the turbulence is fully isotropic: Cx(|m^|) = C||(|m[||) = C(|m'|). Denoting e = |mr|/(C(|m^|)) 
we find: 



OJmR 



1 + 



l-z/Z \2/e-l 



1 



2-e\ 



= COBRi\f 



2-e 



(B5) 



Note that if 6 ^ (there is no accretion thorugh the magnetosphere), a>m ojb, i.e. a solid-body 
rotation without accretion is established (cf. the case one above!). In the case for e = 3/2, a near 
quasi- Keplerian angular rotation law can be realized. We renmind that a similar quasi-Keplerian 
law was obtained in Appendix A above with the use of the Prandtl rule for isotropic turbulent 
viscosity. This was the only solution in that case. Here, in contrast, the quasi-Keplerian law is 
only a partucular case of the general solution obtained with the use of the Wasiutinsky rule for 
anisotropic turbulent viscosity. 

As we have shown in the main text, the case of quasi-Kpeleerian rotation law is not favored 
by observations. So we conclude that most likely near iso-angular-momentum rotational velocity 
distribution with anisotropic turbulence initiated by convection in the shell is realized. We remind 
that in thin accretion discs where the vertical height limits the scale of turbulence, the Prandlt law 
for viscosity works very well (Shakura & Sunyaev, 1973). 
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APPENDIX C: CORRECTIONS TO THE RADIAL TEMPERATURE GRADIENT 

Here we shall estimate how the radial temperature gradient differs from the adiabiatic law due 
to convective motions in the shell. By multiplying Eq. (l25l) by (l/2)(a)„, - a**), we obtain the 
convective heating rate caused by interaction of the shell with the magneto sphere: 

L, = ]^ZMRl{io^ - io*f . (CI) 

Multiplying the same Eq. (|25] ) by of yields the rate of change of the mechanical energy of the 
neutron star 

Lu = ZMRlco*(oj,„ - oj*) . (C2) 
The total energy balance is then 

L, = L, + L,= ^-ZMR\{col - of^) . (C3) 

Note that the obtained formula for L,. is similar to that describing energy release in the boundary 
layer of an accretion disc, see Shakura & Sunyaev (1988). 
The convective energy flux is: 

U ZMRl{co,„-oj*f 

"^' = 4^^ = 8^^^^ • ^"^^^ 

On the other hand, the convective e nergy flux can be related to the entropy gradient (see 



Shakura. Sunyaev & ZilitinkevichI (119781) ) 



qc = -pVcT^, (C5) 
dR 

where S is the specific entropy (per gram). Here Vc is the radial turbulent heat conductivity 

Vc =< uJc >= ChUcR (C6) 

where k ~ R, the convective velocity Uc ~ Cs ~ R~^'^ and C/, is a numerical coefficient of the order 
of one. So 

yc = yciRA)\ — \ ■ (C7) 



Next, we make use of the thermodynamic identity for the specific enthalpy H 
dH \dP. dS 

— = - + T — . (C8) 

dR p dR dR 

We remind the reader that the enthalpy can be written as 

dH = CpdT . 

where 
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\dT 



7 



is the specific heat capacity at constant pressure. Expressing T{dS jdR) from Eq. (ICS I) and making 
use of the hydrostatic equation [Eq. (01)] written as 



dPJp 



n GM 



-tf/(y, nil) , 



dR fimCp 

the thermodynamic identity Eq. (IC8h can be rewritten in the form 



dT_ 
dR 



1 



GM 



i?2 



■iA(r, rtit) - 



Zur{Ra) (Ra 



2v,(Ra) \ R 



*\2 



— ]RUa)^-aj*) 



(C9) 



By definition the adiabatic temperature gradient is determined by the first term on the right side, 
(dT/dR)ad = g/Cp. Equation (IC9I ) can be integrated to find the actual temperature dependence: 



T = — 



GM Zur{Ra)^, 



R 



R 

Ra 



Close to the equilibrium (Id)* = 0), we can use Eq. (IA36I) and write 



r = l 



'P L 



GM 



R 



-ifriy, nit) 



Ur(Ra) 

IChUciRA) 



R 

Ra 



(CIO) 



(Cll) 



This solution shows that in the whole region between R^ and Rg, for slowly rotating pulsars 
(i.e., in which co^ ^ (^k(Ra)), the temperature distribution is close to the adiabatic law with a 
temperature gradient close to the adiabatic one [Eq. (H))]: 
y-\GM 



-if/iy, nit) 



(C12) 



r nR 

In fact, we have only taken into account the energy release due to the frequency difference 
near the magneto sphere. In fact, there can be additional sources of energy in the shell (e.g. the heat 
release during magnetic reconnection and turbulence (see Appendix D), etc.). 



APPENDIX D: DYNAMICS OF A STATIONARY SPHERICALLY-SYMMETRIC 
IDEAL GAS FLOW 

In this Appendix, we write down the gas-dynamic equations of the spherically symmetric ideal 
gas flow onto a Newtonian gravitating center. This problem was considered in the classical pa- 
per by Bondi (1952) for an adibatic accretion flow. Adiabatic gas outflows (stellar winds) were 
studied by Parker (1963). We focus on the role of the cooling/heating processes near the Alvenic 
surface. As we discussed in the main text, at low X-ray luminosities the quasi-static shell is capa- 
ble of removing the angular momentum from the rotating magnetosphere via convective motions. 
As the accretion rate exceeds some critical value, strong Compton cooling causes a free-fall gap 
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to appear above the magnetosphere, and the angular momentum cannot be transferred from the 
magnetosphere to the shell any more. 

The equation of motion Eq. (|A16I) in the absence of viscosity reads: 

du \dP, \dP\ 2(P'-P'J GM 

u — = (Dl) 

dR p dR p dR pR R^ 

Here Pg = pc^Jy is the gas pressure, and P' stands for the pressure due to turbulent pulsations, 
which in general are anisotropic: 

P\\=P< uj >= pm^^c] = yPgml (D2) 

P'^ = 2p< u\ >= 2pm\c] = 2yPgm\ (D3) 

(< >=< u^^ > +2 < u]_ > is the turbulent velocity dispersion, and m\ are the parallel and 
perpendicular turbulent Mach numbers squared). 

From the first law of thermodynamics 
dE Pp dp ^dS 

— = -i-^ + T — (D4) 
dR p dR dR 

where the specific internal energy (per gram) is 

E = CvT= ' (D5) 

r(r - 1) 

the heat capacity is 
H 1 

cv = 7 (D6) 

i^m y - 1 

From the second law of thermodynamics, the specific entropy change can be written through the 
specific heat change rate dQ/dt [erg/s/g] as 

dS dQ dQ/dt ^^^^ 
T — = — = . (D7) 

dR dR u 

Using the mass continuity equation 

M = AnR^pu (D8) 



we find 

I dp 2 1 du^ 



pdR R 2m2 dR ■ 

Using the relation c] = yKT, we finally obtain: 



(D9) 



,2 



1 dCs 

= (r- 1) 

ct dR 



2 1 du 



dQ/dt , 

UCyl 



R 2m2 dR 

Note that this equation can also be derived directly from the ideal gas equation of state written in 
the form 
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where ^ is a constant. 

Using Eq. (|D10I) . the gas pressure gradient can be rewritten in the form: 



(Dll) 



1 dPg 



cj. dQ/dt 2 



CyU 



2 
R 



1 du^ 
2m2 dR 



(D12) 



Plugging Eq. (|D12I) into the equation of motion finally yields: 



1 1 du^ 

iI^Ir 



(ml - ml) GM 



/[u^-cla+yml)] . (D13) 



R 7?2 

Note also that in the strongly anisotropic case where mj = m^ » m]^, the role of turbulence 
increases in comparison with the isotropic case where mj = m]_ = {\/3)mj. 

We can also introduce the Mach number in the flow M = u/cs. Then from Eq. (ID 101) and 
Eq. (ID 131 ) we derive the equation for the Mach number: 



M- dR 

2[(y-l)M^-(y+l)(m^-mi)] _ [M^+y(l+ym^)] _ (y^i)GM \ 

cvT dR R\j j 



(D14) 



R CvT dR 

where we have substituted (dQ/dt) = u{dQ/dR). Equations (IDIOI) . (|D13I) and (|D14|) can be used 
to solve the dynamics of the accretion flow for pairs of independent variables (u, c^), (u, M) or 
(Cv, M). Here, however, we shall only consider the behaviour of the flux near the singular point. 
To this end, we can use Eq. (|D13I) . 

Eq. (ID 131 ) has a singular saddle point where the denominator vanishes: 

M^ = c^(l+rmJ) (D15) 

So must the numerator, from which we find the quadratic equation for the velocity at the singular 
point: 

+ (y - l)mf. + ml^ 



,2 
u — 
R 



1 + 



(dQ/dt\ GM ^ 

- u — = . 

I CvT I i?2 



(D16) 



Remember that in the adiabatic case (dQ/dt = 0) without turbulence at the saddle point we have 
simply 

GM 



2 2 
U = C = 



2R 



(D17) 



We stress that the presence of turbulence increases the velocity at the singular point. For example, 
for 7 = 5/3 we find for strong anisotropic turbulence = c^(l + (5/3)my); for the isotropic 
turbulence the correction is smaller: = cj(l + (5/9)m^). The transition through the sound speed 
(the sonic point where = c^) lies above the saddle point due to turbulence, and there is no 
singularity in the sonic point. 
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First let us determine the turbulence heating rate in the quasi-static shell (dQ/dt)^: 



dQ 
dt 



l<u, > 



2 tt 

where the characteristic time of the turbulent heating is 
R R 



tt = at— = at 

Ut nitc. 



(D18) 



(D19) 



with at being a dimensionless constant characterizing the turbulent dissipation energy rate and the 
turbulent Mach number is mj = m?: + 2m\. The turbulent heating rate can thus be written as 



dQ 
dt 



cl 



-m: 



latR "* ■ 

In the case of Compton cooling we have 



dQ 
dt 



cyjT - r,) 

tc 



(D20) 



(D21) 



where tc is the Compton cooling time [Eq. (I34l)1. 
Eq. (ID 161) can now be written in the form: 
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(D22) 



u 2a,R ytc R^ 

As our problem is one of accretion, the sign of the velocity u = dR/dt is negative, so below 
we shall write u = -\u\. Then for the absolute value of the velocity at the singular point where 
Cs/M = -1/(1 + yniy^'^ we have the quadratic equation: 



22 
u — 
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2 \ 



'1 + (y- \)ml + m]_ 



1 + ym^ 
ase, th( 

Ra - TjT) 



+ u 



1 y(y-lK \u\(l-TJT) GM 



(l+ymJ)i/2 2atR 
In this case, the solution to Eq. (ID 161) reads: 

R R^H - TJTf 
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(D23) 



(D24) 



(D25) 



where we have introduced the dimensionless factor 
1 + (y - l)ml + m\ y{y - V)(m\ + 2m^J'''^ 

l+ymf^ 4at(l+ymjy/^ 

In the case of isotropic turbulence where my = = 1/ V3, m, = 1, for y = 5/3 the factor 
A X 1.23, and in the case of strongly anisotropic turbulence where my = l,mx = 0, = 1, this 
factor is A « 0.8. 

In units of the free-fall velocity the solution Eq. (|D24I) reads: 



^ \u\ (1-TJT) 

f(u) = — = — 

Uff AyA 



1 ^ (1 - TJTf 
A ^ Ay^A^ 



'21 

tc 



1/2 



(D26) 



With Compton cooling, the temperature changes exponentially: 
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T = T, + {T,r - T,)e-'l'^ (D27) 

(see the main text). When cooUng is slow, tffjtc 1, the critical point lies inside the Alfven 
surface, i.e. no transition through the critical point occurs in the flow before it meets the magneto- 
sphere, and in this case we expect that the accretion settling regime can be realized. If this point 
lies abobe the Alfven surface, the velocity of the flow can become supersonic above the magneto- 
sphere, and one can expect the formation of a shock. Both turbulence and rapid cooling shifts the 
location of the critical point upward in the flow. 

In the case of rapid cooling, tff/tc » 1, T ^ T^, so again m/m// » 1/2, but the critical 
point rises above the Alfven surface, so a free-fall gap above the magnetosphere appears. The 
ratio f{u) = \u\/uff reaches a maximum at tff/tc ~ 0.46 (assuming a typical ratio Tcr/T^ = 10), 
and depending on the value of A = 0.8 ^ 1.23 (anisotropic or isotropic turbulence) it equals 
/(m) = 0.5 - 0.6. 

ACKNOWLEDGMENTS 

We are very grateful to the anonymous referee for careful reading of the manuscript and many 
useful comments. The authors also acknowledge Dr. V. Suleimanov (lAAT) for discussions and 
valuable notes, Dr. V. Doroshenko (lAAT) for courtesly providing the torque-luminosity plots for 
GX 301-2 and Vela X-1, and Ms. A. Gonzalez- Galan for providing data on GX 1-1-4 prior to 
publication. NIS thanks the Max-Planck Institute for Astrophysics (Garching) for the hospitality. 
The work by NIS, KAP and AYK is supported by RFBR grants 09-02-00032, 12-02-00186, and 
10-02-00599. LH is supported by a grant from the Wenner-Gren foundations. 

REFERENCES 

Arons J., & Lea S. M. 1976a ApJ, 207, 914 

Arons J., & Lea S. M. 1976b ApJ, 210, 792 

Bildsten L., et al., 1997, ApJS, 113, 367 

Bisnovatyi-Kogan G. S. 1991, A&A, 245, 528 

Bochkarev N. G., Karitskaja E. A., Shakura N. I. 1975, SvAL 1, 237 

Bondi H. 1952, MNRAS, 112, 195 

Burnard D.J., Arons J., Lea S. M. 1983, ApJ, 266, 175 

Chakrabarty D., et al. 1997, ApJL, 481, LlOl 



Quasi-spherical accretion 47 

Chashkina A., & Popov S.B., 2011, submitted to New Astronomy 
Cowie L. L., McKee C. R, & Ostriker J. P. 1981, ApJ, 247, 908 
Davidsen A., Malina R., & Bowyer S. 1977, ApJ, 21 1, 866 
Davies R. E., & Pringle J. E. 1981, MNRAS, 196, 209 
de Kool M., & Anzer U. 1993, MNRAS, 262, 726 

Doroshenko V., Santangelo A., Suleimanov V., Kreykenbohm I., Staubeit R., Ferrigno C, & 

Klochkov D. 2010, A&A, 515, AlO 
Dotani,T., Kii T., Nagase R, Makishima K., Ohashi T., Sakao T., Koyama K., & Tuohy, 1. R. 

1989, PASJ, 41, 427 
Ducci L., Sidoli L. & Paizis A. 2010, MNRAS, 408, 1540 
Eisner R. R, & Lamb R K. 1977, ApJ, 215, 897 

Ringer M. et al., http://gammaray.nsstc.nasa.gov/gbm/science/pulsars/lightcurves/gxlp4.html 
Rryxell B. A., & Taam R. E. 1988, ApJ, 335, 862 
Ghosh P, & Lamb R K. 1979, ApJ, 234, 296 
Gonzalez-Galan et al., submitted to AA llarXiv: 1105. 190711 
Hatchett S., BufF J., & McCray R. 1976, ApJ, 206, 847 

Hinkle K. H., Fekel R C, Joyce R. R., Wood P R., Smith V. V., Lebzelter T. 2006, ApJ, 641, 479 
Ho C, Taam R. E., Fryxell B. A. et al. 1989, MNRAS, 238, 1447 
HuntR., 1971, MNRAS, 154, 141 

lUarionov A. R, & Kompaneets D. A. 1990, MNRAS, 247, 219 

lUarionov A. R& Sunyaev R. A., 1975, A&A, 39, 185 

Kaper L., van der Meer A., & Najarro R 2006, A&A, 457, 595 

Kluzniak W., & Rappaport S. 2007, ApJ, 671, 1990 

Koh D. T., et al. 1997, ApJ, 479, 933 

Kompaneets A.S, 1956, ZhETP, 31, 876 

Kreykenbohm I., Wilms J., Cobum W., Kuster M., Rothschild R. E., Heindl W. A., Kretschmar 

P, Staubert R. 2004, A&A, 427, 975 
La Barbera A., Segreto A., Santangelo A., Kreykenbohm I., Orlandini M. 2005, A&A, 438, 617 
Landau D., Lifshitz E. M., Hydrodynamics, Nauka, Moscow, 1986 
Lovelace R. V. E., Romanova M. M., Bisnovatyi-Kogan G. S. 1995, MNRAS, 275, 244 
Marykutty J., Biswajit P, Jincy D., Kavila I. 2010, MNRAS, 407, 285 
Makishima K., et al. 1988, Nature, 333, 746 

Nagase, R, Hayakawa, S., Sato, N., Masai, K., & Inoue, H. 1986, PASJ, 38, 547 



48 A^. Shakura et al. 
NagaseF. 1989, PAS J, 41, 1 

Nelson R.W., Bildsteen L., Chakrabarty D., et al. 1997, ApJ, 488, LI 17 

Parker, E. 1963, Interplanetary dynamical processes (New York, Interscience Publishers) 

Pravdo S. H., & Ghosh R 2001, ApJ, 554, 383 

Pringle J.E., & Rees M.J. 1972, A&A, 21, 1 

Quaintrell, H., Norton, A. J., Ash, T. D. C, Roche, P, Willems, B., Bedding, T. R., Baldry, 1. K., 

& Fender, R. P 2003, A&A, 401, 313 
Raymond J. C, Cox D. P, & Smith B. W. 1976, ApJ, 204, 290 
Rappaport, S. 1975, lAU Circ, 2869, 2 
RufFert M. 1997, A&A, 317, 793 
RufFert M. 1999, A&A, 346, 861 
Shakura N.I., & Sunyaev R.A., 1973, A&A, 24, 337 
Shakura N.L, Sunyaev R.A., & Zilitinkevich 1978, A&A, 62, 179 
Shakura N.L, & Sunyaev R.A. 1988, Adv. Space Res., 8, 135 

Staubert R. 2003, in Chinese Journal of Astronomy and Astrophysics Supplement, 3, 270 
Sunyaev R.A. 1978, in: Physics and astrophysics of neutron stars and black holes, Amsterdam, 
North Holland PubUshing Co., p. 697 
Sunyaev R.A., & Shakura N.L, 1977, SvAL, 3, 138 
Tarter C. B., Tucker W. H., & Salpeter E. E. 1969, ApJ, 156, 943 
van Kerkwijk M. H., van Paradijs J., Zuiderwijk E. J., et al. 1995, A&A, 303, 483 
Wasiutinski J., 1946, Studies in Hydrodynamics and Structure of Stars and Planets, Oslo. 
Watanabe, S., et al. 2006, ApJ, 651, 421 

White N. E., Mason K. O., Huckle H. E., Charles P A., Sanford P W. 1976, ApJL, 209, LI 19 
Weymann R. 1965, Phys. Fluids, 8, 2112 



